# Directory and Data Import

setwd("/home/joshk/git_repositories/BrookTrout_TG/Data")
Trout = read.csv("BrookTrout.csv")
    Trout = subset(Trout, !is.na(Mother.ID))
    Trout = subset(Trout, !is.na(Father.ID))
Trout$Offspring.Acc = as.factor(Trout$Offspring.Acc)
Trout$Mother.Acc = as.factor(Trout$Mother.Acc)
Trout$Father.Acc = as.factor(Trout$Father.Acc)
Trout$Mother.ID = as.factor(Trout$Mother.ID)
Trout$Father.ID = as.factor(Trout$Father.ID)

options(na.action = "na.fail")

## Packages

library('easypackages')
libraries('car','dotwhisker',"dplyr","ggplot2","grid","gridExtra",
          "knitr","lme4","magrittr", "MuMIn", "nlme","purrr","rmarkdown")

Sectioning Data

Mass.Spec = Trout[,c(1:6, 11:17)]
N.Mass.Spec = Trout[,c(1:10, 14:17)]

Assigning ggplot2 themes for figures

my.theme = theme(panel.grid.minor = element_blank(), 
axis.title = element_text(size = 14, family = "Noto Sans"), 
axis.text = element_text(size = 14, colour = "black", family = "Noto Sans"), 
axis.title.y = element_text(angle = 0, vjust = 0.5), panel.grid.major = 
element_line(colour = "grey75"), legend.title = element_text(size = 14, 
colour = "black", family = "Noto Sans"), legend.text = element_text(size = 12, 
colour = "black", family = "Noto Sans"))

my.theme.dw = theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), 
axis.title = element_text(size = 12, family = "Noto Sans"), 
axis.text = element_text(size = 12, colour = "black", family = "Noto Sans"), 
legend.title = element_text(size = 12, 
colour = "black", family = "Noto Sans"), legend.text = element_text(size = 12, 
colour = "black", family = "Noto Sans"))

my.theme.diag = theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), 
axis.title = element_text(size = 8, family = "Noto Sans"), 
axis.text = element_text(size = 8, colour = "black", family = "Noto Sans"), 
legend.title = element_text(size = 8, 
colour = "black", family = "Noto Sans"), legend.text = element_text(size = 8, 
colour = "black", family = "Noto Sans"))

Assessing distribution of response variables

{
  RMR.MS.Hist = ggplot(data = Mass.Spec, aes(x = RMR.MS)) + 
  geom_histogram(colour = "black", fill = "royalblue") + my.theme + theme_bw() +
  geom_vline(xintercept=with(na.omit(Mass.Spec), mean(RMR.MS)), size = 1.25, 
  colour="darkred", linetype = "dotted") + xlab("Resting Metabolic \n Rate (mg O2/h)")

  MMR.MS.Hist = ggplot(data = Mass.Spec, aes(x = MMR.MS)) + 
  geom_histogram(colour = "black", fill = "royalblue") + my.theme + theme_bw() +
  geom_vline(xintercept=with(na.omit(Mass.Spec), mean(MMR.MS)), size = 1.25, 
  colour="darkred", linetype = "dotted") + xlab("Maximum Metabolic \n Rate (mg O2/h)")

  AAS.MS.Hist = ggplot(data = Mass.Spec, aes(x = AAS.MS)) + 
  geom_histogram(colour = "black", fill = "royalblue") + my.theme + theme_bw() +
  geom_vline(xintercept=with(na.omit(Mass.Spec), mean(AAS.MS)), size = 1.25, 
  colour="darkred", linetype = "dotted") + xlab("Absolute Aerobic \n Scope")

  RMR.Hist = ggplot(data = N.Mass.Spec, aes(x = RMR)) + 
  geom_histogram(colour = "black", fill = "royalblue") + my.theme + theme_bw() +
  geom_vline(xintercept=with(na.omit(N.Mass.Spec), mean(RMR)), size = 1.25, 
  colour="darkred", linetype = "dotted") + xlab("Resting Metabolic \n Rate (mg O2/h)")

  MMR.Hist = ggplot(data = N.Mass.Spec, aes(x = MMR)) + 
  geom_histogram(colour = "black", fill = "royalblue") + my.theme + theme_bw() +
  geom_vline(xintercept=with(na.omit(N.Mass.Spec), mean(MMR)), size = 1.25, 
  colour="darkred", linetype = "dotted") + xlab("Maximum Metabolic \n Rate (mg O2/h)")

  AAS.Hist = ggplot(data = N.Mass.Spec, aes(x = AAS)) + 
  geom_histogram(colour = "black", fill = "royalblue") + my.theme + theme_bw() +
  geom_vline(xintercept=with(na.omit(N.Mass.Spec), mean(AAS)), size = 1.25, 
  colour="darkred", linetype = "dotted") + xlab("Absolute Aerobic \n Scope")
}

grid.arrange(RMR.MS.Hist, MMR.MS.Hist, AAS.MS.Hist, 
RMR.Hist, MMR.Hist, AAS.Hist, nrow = 2, 
top = "Row 1 : Mass Specific, Row 2: Non-Mass Specific")  

Mass.hist = ggplot(data = Trout, aes(x = Mass)) + 
geom_histogram(colour = "black", fill = "royalblue") + my.theme + theme_bw() +
geom_vline(xintercept=with(Trout, mean(Mass)), size = 1.25, 
colour="darkred", linetype = "dotted") + xlab("Mass")

CTM.hist = ggplot(data = Trout, aes(x = CTM)) + 
geom_histogram(colour = "black", fill = "royalblue") + my.theme + theme_bw() +
geom_vline(xintercept=with(subset(N.Mass.Spec, !is.na(CTM)), mean(CTM)), size = 1.25, 
colour="darkred", linetype = "dotted") + xlab("Critical Thermal Maximum (Degrees C)")

grid.arrange(Mass.hist, CTM.hist, nrow = 1)

## Mass appears to fit a Gaussian distibution, but mass specific
## parameters do not. Critical thermal maximum not awfully skewed, and all 
## non-mass specific parameters meet gaussian distributions.

Scaling with Mass

RMR.by.Mass = ggplot(data = N.Mass.Spec, aes(x = Mass, y = RMR)) + 
    my.theme + theme_bw() + geom_point(size = 3, colour = "cadetblue", alpha = 0.3) + xlab("Mass (g)") + 
    geom_smooth(method = lm, colour = "black") + ylab("Resting Metabolic Rate \n (mg O2/h)")

MMR.by.Mass = ggplot(data = N.Mass.Spec, aes(x = Mass, y = AAS)) + 
    my.theme + theme_bw() + geom_point(size = 3, colour = "palegreen3", alpha = 0.3) + xlab("Mass (g)") +
    geom_smooth(method = lm, colour = "black") + ylab("Maximum Metabolic Rate \n (mg O2/h)")    

AAS.by.Mass = ggplot(data = N.Mass.Spec, aes(x = Mass, y = RMR)) + 
    my.theme + theme_bw() + geom_point(size = 3, colour = "thistle", alpha = 0.3) + xlab("Mass (g)") +
    geom_smooth(method = lm, colour = "black") + ylab("Absolute Aerobic Scope") 

CTM.by.Mass = ggplot(data = N.Mass.Spec, aes(x = Mass, y = CTM)) + 
    my.theme + theme_bw() + geom_point(size = 3, colour = "gold2", alpha = 0.3) + xlab("Mass (g)") +
    geom_smooth(method = lm, colour = "black") + ylab("Critical Thermal Maximum \n (Degrees C)")    

grid.arrange(RMR.by.Mass, MMR.by.Mass, AAS.by.Mass, CTM.by.Mass, 
nrow = 2, top = "Metabolic Rate and Scope by Mass")  

# No clear outliers, but critical thermal maximum and resting metabolic rate heterogenious
# by mass. Weighting will be necessary if mass-specific responses not used.

Assessing Outliers with Dotplots

RMR.NMS.Data = subset(N.Mass.Spec, !is.na(RMR))
RMR.dotplot = ggplot(data = RMR.NMS.Data, aes(x = RMR, y = 
    (1:nrow(RMR.NMS.Data)))) + my.theme + theme_bw() + geom_point(size = 
    3, colour = "gold2", alpha = 0.5) + xlab("Resting Metabolic Rate (mg O2/h)") +
    ylab("Row Number") + annotate("rect", xmin = (with(RMR.NMS.Data, mean(RMR) - 3*sd(RMR))),
    xmax = (with(RMR.NMS.Data, mean(RMR) + 3*sd(RMR))), ymin = 0, ymax = nrow(RMR.NMS.Data),
    alpha = .2)
MMR.NMS.Data = subset(N.Mass.Spec, !is.na(MMR))
MMR.dotplot = ggplot(data = MMR.NMS.Data, aes(x = MMR, y = 
    (1:nrow(MMR.NMS.Data)))) + my.theme + theme_bw() + geom_point(size = 
    3, colour = "royalblue", alpha = 0.5) + xlab("Maximum Metabolic Rate (mg O2/h)") +
    ylab("Row Number") + annotate("rect", xmin = (with(MMR.NMS.Data, mean(MMR) - 3*sd(MMR))),
    xmax = (with(MMR.NMS.Data, mean(MMR) + 3*sd(MMR))), ymin = 0, ymax = nrow(MMR.NMS.Data),
    alpha = .2)
AAS.NMS.Data = subset(N.Mass.Spec, !is.na(AAS))
AAS.dotplot = ggplot(data = AAS.NMS.Data, aes(x = AAS, y = 
    (1:nrow(AAS.NMS.Data)))) + my.theme + theme_bw() + geom_point(size = 
    3, colour = "darkseagreen2", alpha = 0.5) + xlab("Absolute Aerobic Scope") +
    ylab("Row Number") + annotate("rect", xmin = (with(AAS.NMS.Data, mean(AAS) - 3*sd(AAS))),
    xmax = (with(AAS.NMS.Data, mean(AAS) + 3*sd(AAS))), ymin = 0, ymax = nrow(AAS.NMS.Data),
    alpha = .2)
CTM.NMS.Data = subset(N.Mass.Spec, !is.na(CTM))
CTM.dotplot = ggplot(data = CTM.NMS.Data, aes(x = CTM, y = 
    (1:nrow(CTM.NMS.Data)))) + my.theme + theme_bw() + geom_point(size = 
    3, colour = "black", alpha = 0.5) + xlab("Critical Thermal Maximum (Degrees C)") +
    ylab("Row Number") + annotate("rect", xmin = (with(CTM.NMS.Data, mean(CTM) - 3*sd(CTM))),
    xmax = (with(CTM.NMS.Data, mean(CTM) + 3*sd(CTM))), ymin = 0, ymax = nrow(CTM.NMS.Data),
    alpha = .2)

grid.arrange(RMR.dotplot, MMR.dotplot, AAS.dotplot, CTM.dotplot, 
nrow = 2, top = "Non-Mass Specific Response Dotplots: \n Grey boxed represent mean +/-  
    three standard deviations")  

## Outliers not extremely distinct, so retained but monitored below.

## Repeating dotplots for mass-specific data.

RMR.MS.Data = subset(Mass.Spec, !is.na(RMR.MS))
RMR.dotplot = ggplot(data = RMR.MS.Data, aes(x = RMR.MS, y = 
    (1:nrow(RMR.MS.Data)))) + my.theme + theme_bw() + geom_point(size = 
    3, colour = "gold2", alpha = 0.5) + xlab("Resting Metabolic Rate (mg O2/kg/h)") +
    ylab("Row Number") + annotate("rect", xmin = (with(RMR.MS.Data, mean(RMR.MS) - 3*sd(RMR.MS))),
    xmax = (with(RMR.MS.Data, mean(RMR.MS) + 3*sd(RMR.MS))), ymin = 0, ymax = nrow(RMR.MS.Data),
    alpha = .2)
MMR.MS.Data = subset(Mass.Spec, !is.na(MMR.MS))
MMR.dotplot = ggplot(data = MMR.MS.Data, aes(x = MMR.MS, y = 
    (1:nrow(MMR.MS.Data)))) + my.theme + theme_bw() + geom_point(size = 
    3, colour = "royalblue", alpha = 0.5) + xlab("Maximum Metabolic Rate (mg O2/kg/h)") +
    ylab("Row Number") + annotate("rect", xmin = (with(MMR.MS.Data, mean(MMR.MS) - 3*sd(MMR.MS))),
    xmax = (with(MMR.MS.Data, mean(MMR.MS) + 3*sd(MMR.MS))), ymin = 0, ymax = nrow(MMR.MS.Data),
    alpha = .2)
AAS.MS.Data = subset(Mass.Spec, !is.na(AAS.MS))
AAS.dotplot = ggplot(data = AAS.MS.Data, aes(x = AAS.MS, y = 
    (1:nrow(AAS.MS.Data)))) + my.theme + theme_bw() + geom_point(size = 
    3, colour = "darkseagreen2", alpha = 0.5) + xlab("Absolute Aerobic Scope") +
    ylab("Row Number") + annotate("rect", xmin = (with(AAS.MS.Data, mean(AAS.MS) - 3*sd(AAS.MS))),
    xmax = (with(AAS.MS.Data, mean(AAS.MS) + 3*sd(AAS.MS))), ymin = 0, ymax = nrow(AAS.MS.Data),
    alpha = .2)
Details = ggplot(data = data.frame("x" = c(1:100), "y" = c(1:100)), aes(x = x, y = y)) +
theme_void() + theme(legend.position="none") + annotate("text", x = 25, y = 70, label = 
"Grey boxes represent mean +/-  \n three standard deviations")

    grid.arrange(RMR.dotplot, MMR.dotplot, AAS.dotplot, Details,  
    nrow = 2, top = "Mass Specific Response Dotplots")

# One outlier, though not visually concerning. Retaining and monitoring below.

##Non-Mass Specific Linear Models

####Note that sires are not unique to dams; therefore random effects should be crossed

RMR.model = lmer(RMR ~ Mother.Acc*Father.Acc*Offspring.Acc + Mass + (1|Mother.ID) +(1|Father.ID), 
    REML = FALSE, data = RMR.NMS.Data)  
MMR.model = lmer(MMR ~ Mother.Acc*Father.Acc*Offspring.Acc + Mass + (1|Mother.ID) +(1|Father.ID),
    REML = FALSE, data = MMR.NMS.Data)
AAS.model = lmer(AAS ~ Mother.Acc*Father.Acc*Offspring.Acc + Mass + (1|Mother.ID) +(1|Father.ID), 
    REML = FALSE, data = AAS.NMS.Data)
CTM.model = lmer(CTM ~ Mother.Acc*Father.Acc*Offspring.Acc + Mass + (1|Mother.ID) +(1|Father.ID), 
    REML = FALSE, weights = Mass, data = CTM.NMS.Data)

# Sample sizes

nrow(RMR.NMS.Data)
## [1] 186
nrow(MMR.NMS.Data)
## [1] 224
nrow(AAS.NMS.Data)
## [1] 170
nrow(CTM.NMS.Data)
## [1] 216

###Quick Assessment of Dispersion

{
  RMR.Dev = abs(summary(RMR.model)$AICtab[[4]]/summary(RMR.model)$AICtab[[5]])
  MMR.Dev = abs(summary(MMR.model)$AICtab[[4]]/summary(MMR.model)$AICtab[[5]])
  AAS.Dev = abs(summary(AAS.model)$AICtab[[4]]/summary(AAS.model)$AICtab[[5]])
  CTM.Dev = abs(summary(CTM.model)$AICtab[[4]]/summary(CTM.model)$AICtab[[5]])
}
print(paste("RMR:", RMR.Dev, "MMR:", MMR.Dev, "AAS:", AAS.Dev, "CTM:", CTM.Dev, sep = " "))
## [1] "RMR: 0.363211108446256 MMR: 0.0141738136619573 AAS: 0.271903624981455 CTM: 0.0593597551129788"
# All highly over-dispersed

###Repeating with log Transformed Response Variable and Mass; Save for AAS.

RMR.model = lmer(log(RMR) ~ Mother.Acc*Father.Acc*Offspring.Acc + log(Mass) + (1|Mother.ID) +(1|Father.ID), 
    REML = FALSE, data = RMR.NMS.Data)  
MMR.model = lmer(log(MMR) ~ Mother.Acc*Father.Acc*Offspring.Acc + log(Mass) + (1|Mother.ID) +(1|Father.ID),
    REML = FALSE, data = MMR.NMS.Data)
AAS.model = lmer(AAS ~ Mother.Acc*Father.Acc*Offspring.Acc + log(Mass) + (1|Mother.ID) +(1|Father.ID), 
    REML = FALSE, data = AAS.NMS.Data)
CTM.model = lmer(CTM ~ Mother.Acc*Father.Acc*Offspring.Acc + log(Mass) + (1|Mother.ID) +(1|Father.ID), 
    REML = FALSE, weights = Mass, data = CTM.NMS.Data)

# Double-checking dispersion 
{
  RMR.Dev = abs(summary(RMR.model)$AICtab[[4]]/summary(RMR.model)$AICtab[[5]])
  MMR.Dev = abs(summary(MMR.model)$AICtab[[4]]/summary(MMR.model)$AICtab[[5]])
  AAS.Dev = abs(summary(AAS.model)$AICtab[[4]]/summary(AAS.model)$AICtab[[5]])
  CTM.Dev = abs(summary(CTM.model)$AICtab[[4]]/summary(CTM.model)$AICtab[[5]])
}
print(paste("RMR:", RMR.Dev, "MMR:", MMR.Dev, "AAS:", AAS.Dev, "CTM:", CTM.Dev, sep = " "))
## [1] "RMR: 1.29944000445268 MMR: 0.0456298181606057 AAS: 0.285021278373539 CTM: 0.0595168880242838"
## Slightly better.

# Checking sample size for each grouping, and calculating power

RMR.NMS.Data %>% group_by(Offspring.Acc, Father.Acc, Mother.Acc) %>% 
  summarize(N = n()) %>% as.data.frame() %>% summarize(Mean = mean(N), SD = sd(N))
##    Mean       SD
## 1 23.25 5.700877
require('pwr')

num_df = nrow(summary(RMR.model)$coefficients) - 1
den_df = nrow(RMR.NMS.Data) - (num_df + 1)
pwr.f2.test(u = num_df, v = den_df, f2 = 0.05, sig.level = 0.05)
## 
##      Multiple regression power calculation 
## 
##               u = 8
##               v = 177
##              f2 = 0.05
##       sig.level = 0.05
##           power = 0.5284571

##Conducting Model Iterations and Sorting by AIC

AIC.Results.NMS = vector('list', 4)
models = c("RMR.model", "MMR.model", "AAS.model", "CTM.model")

# Looping through models, running with all combinations of predictors,
# then sorting by AICc

for (i in 1:length(models)){
  AIC.Results.NMS[[i]] = subset(dredge(get(models[i]), rank = "AICc", 
  extra = list(R2 = function(x) r.squaredGLMM(x))), delta < 5)
}

Select.Models.NMS = c()
Delta.NMS = c()
Response = c()
Marg.R2 = c() 
Cond.R2 = c()

# Looping through models and collecting metrics of fit (i.e. marginal and conditional r-squared values)

for (j in 1:length(AIC.Results.NMS)){
  for (i in 1:length(attr(AIC.Results.NMS[[j]], "model.calls"))){
    Select.Models.NMS[i] = gsub(".*~ ", "", paste(attr(AIC.Results.NMS[[j]], "model.calls")[[i]])[2])
    Delta.NMS[i] = AIC.Results.NMS[[j]]$delta[i]
    Response[i] = gsub(" ~.*", "", paste(attr(AIC.Results.NMS[[j]], "model.calls")[[i]])[2])
    Marg.R2[i] = AIC.Results.NMS[[j]]$R21[i]
    Cond.R2[i] = AIC.Results.NMS[[j]]$R22[i]
    
    if (i == length(attr(AIC.Results.NMS[[j]], "model.calls"))){
      Table = data.frame("Reponse.Type" = "Mass Independant", "Response" = Response, 
    "Delta.AIC" = Delta.NMS, "Marginal.R2" = Marg.R2, "Conditional.R2" = Cond.R2,
    "Equation" = Select.Models.NMS)
    assign(paste("AIC.Selection", "NMS", j, sep = "_"), Table)
    }
  }
  Select.Models.NMS = c()
  Delta.NMS = c()
  Response = c()
  Marg.R2 = c() 
  Cond.R2 = c()
}

AIC.Selection.NMS = rbind(AIC.Selection_NMS_1, AIC.Selection_NMS_2, AIC.Selection_NMS_3, AIC.Selection_NMS_4)

##Conducting Mass-Specific Linear Models

RMR.MS.model = lmer(RMR.MS ~ Mother.Acc*Father.Acc*Offspring.Acc + (1|Mother.ID) + (1|Father.ID), 
            REML = FALSE, data = RMR.MS.Data)
MMR.MS.model = lmer(MMR.MS ~ Mother.Acc*Father.Acc*Offspring.Acc + (1|Mother.ID) + (1|Father.ID), 
            REML = FALSE, data = MMR.MS.Data)
AAS.MS.model = lmer(AAS.MS ~ Mother.Acc*Father.Acc*Offspring.Acc + (1|Mother.ID) + (1|Father.ID), 
            REML = FALSE, data = AAS.MS.Data)

###Quickly Assessing Dispersion.

{
  RMR.MS.Dev = abs(summary(RMR.MS.model)$AICtab[[4]]/summary(RMR.MS.model)$AICtab[[5]])
  MMR.MS.Dev = abs(summary(MMR.MS.model)$AICtab[[4]]/summary(MMR.MS.model)$AICtab[[5]])
  AAS.MS.Dev = abs(summary(AAS.MS.model)$AICtab[[4]]/summary(AAS.MS.model)$AICtab[[5]])
}

print(paste("RMR:", RMR.MS.Dev, "MMR:", MMR.MS.Dev, "AAS:", AAS.MS.Dev, sep = " "))
## [1] "RMR: 12.8500706198959 MMR: 12.1844411795225 AAS: 12.1235924890615"
# All very over-dispersed. Log-transform response and repeat.

RMR.MS.model = lmer(log(RMR.MS) ~ Mother.Acc*Father.Acc*Offspring.Acc + (1|Mother.ID) + (1|Father.ID), 
            REML = FALSE, data = RMR.MS.Data)
MMR.MS.model = lmer(log(MMR.MS) ~ Mother.Acc*Father.Acc*Offspring.Acc + (1|Mother.ID) + (1|Father.ID), 
            REML = FALSE, data = MMR.MS.Data)
AAS.MS.model = lmer(log(AAS.MS) ~ Mother.Acc*Father.Acc*Offspring.Acc + (1|Mother.ID) + (1|Father.ID), 
            REML = FALSE, data = AAS.MS.Data)

{
  RMR.MS.Dev = abs(summary(RMR.MS.model)$AICtab[[4]]/summary(RMR.MS.model)$AICtab[[5]])
  MMR.MS.Dev = abs(summary(MMR.MS.model)$AICtab[[4]]/summary(MMR.MS.model)$AICtab[[5]])
  AAS.MS.Dev = abs(summary(AAS.MS.model)$AICtab[[4]]/summary(AAS.MS.model)$AICtab[[5]])
}

print(paste("RMR:", RMR.MS.Dev, "MMR:", MMR.MS.Dev, "AAS:", AAS.MS.Dev, sep = " "))
## [1] "RMR: 1.20661839775676 MMR: 0.372623478173389 AAS: 0.407045832769152"
# Significantly better.

##Again, Conducting Model Iterations and Sorting by AIC.

AIC.Results.MS = vector('list', 3)
models.MS = c("RMR.MS.model", "MMR.MS.model", "AAS.MS.model")

# Again, looping through models, re-running with all combinations of predictors, then calculating AICc

for (i in 1:length(models.MS)){
  AIC.Results.MS[[i]] = subset(dredge(get(models.MS[i]), rank = "AICc", 
  extra = list(R2 = function(x) r.squaredGLMM(x))), delta < 5)
}

Select.Models.MS = c()
Delta.MS = c()
Response.MS = c()
Marg.R2.MS = c() 
Cond.R2.MS = c()

# Looping through models and pulling out metrics of fit

for (j in 1:length(AIC.Results.MS)){
  for (i in 1:length(attr(AIC.Results.MS[[j]], "model.calls"))){
    Select.Models.MS[i] = gsub(".*~ ", "", paste(attr(AIC.Results.MS[[j]], "model.calls")[[i]])[2])
    Delta.MS[i] = AIC.Results.MS[[j]]$delta[i]
    Response.MS[i] = gsub(" ~.*", "", paste(attr(AIC.Results.MS[[j]], "model.calls")[[i]])[2])
    Marg.R2.MS[i] = AIC.Results.MS[[j]]$R21[i]
    Cond.R2.MS[i] = AIC.Results.MS[[j]]$R22[i]
    if (i == length(attr(AIC.Results.MS[[j]], "model.calls"))){
      Table = data.frame("Reponse.Type" = "Mass-Specific", "Response" = Response.MS, 
      "Delta.AIC" = Delta.MS, "Marginal.R2" = Marg.R2.MS, "Conditional.R2" = Cond.R2.MS,
      "Equation" = Select.Models.MS)
      assign(paste("AIC.Selection", "MS", j, sep = "_"), Table)
    }
  }
  Select.Models.MS = c()
  Delta.MS = c()
  Response.MS = c()
  Marg.R2.MS = c() 
  Cond.R2.MS = c()
}

AIC.Selection.MS = rbind(AIC.Selection_MS_1, AIC.Selection_MS_2, AIC.Selection_MS_3)

# Stacking all top models

AIC.Selection = rbind(AIC.Selection.MS, AIC.Selection.NMS)
print(AIC.Selection, right = F)
##    Reponse.Type     Response    Delta.AIC Marginal.R2 Conditional.R2
## 1  Mass-Specific    log(RMR.MS) 0.0000000 0.09689341  0.23258538    
## 2  Mass-Specific    log(RMR.MS) 0.4320103 0.05317331  0.24043771    
## 3  Mass-Specific    log(RMR.MS) 1.2633111 0.11949983  0.22942058    
## 4  Mass-Specific    log(RMR.MS) 1.6197476 0.07652765  0.23332101    
## 5  Mass-Specific    log(RMR.MS) 2.1592505 0.09688383  0.23257954    
## 6  Mass-Specific    log(RMR.MS) 3.4158810 0.11948859  0.22852419    
## 7  Mass-Specific    log(RMR.MS) 3.4339163 0.11966268  0.23078182    
## 8  Mass-Specific    log(RMR.MS) 3.4465742 0.11949854  0.22941555    
## 9  Mass-Specific    log(RMR.MS) 3.7625995 0.07663967  0.23491827    
## 10 Mass-Specific    log(MMR.MS) 0.0000000 0.12081940  0.14201563    
## 11 Mass-Specific    log(MMR.MS) 1.1337870 0.12482063  0.14707607    
## 12 Mass-Specific    log(MMR.MS) 1.6198084 0.07304761  0.14942805    
## 13 Mass-Specific    log(MMR.MS) 1.6632302 0.12283706  0.14581838    
## 14 Mass-Specific    log(MMR.MS) 2.7260488 0.07701051  0.15563521    
## 15 Mass-Specific    log(MMR.MS) 2.7981663 0.12640953  0.14714116    
## 16 Mass-Specific    log(MMR.MS) 2.8091275 0.12687437  0.15094298    
## 17 Mass-Specific    log(MMR.MS) 3.1451250 0.12529971  0.14712299    
## 18 Mass-Specific    log(MMR.MS) 4.4679396 0.07829618  0.15562601    
## 19 Mass-Specific    log(MMR.MS) 4.5445400 0.12823397  0.15073801    
## 20 Mass-Specific    log(MMR.MS) 4.7897724 0.12703622  0.14722626    
## 21 Mass-Specific    log(MMR.MS) 4.8318359 0.12737623  0.15100909    
## 22 Mass-Specific    log(AAS.MS) 0.0000000 0.06030224  0.06030224    
## 23 Mass-Specific    log(AAS.MS) 1.6058206 0.06331924  0.06331924    
## 24 Mass-Specific    log(AAS.MS) 1.6547206 0.06304827  0.06304827    
## 25 Mass-Specific    log(AAS.MS) 2.7150419 0.06921105  0.06921105    
## 26 Mass-Specific    log(AAS.MS) 3.2960644 0.06600676  0.06600676    
## 27 Mass-Specific    log(AAS.MS) 3.8171171 0.06312377  0.06312377    
## 28 Mass-Specific    log(AAS.MS) 4.4085037 0.07201229  0.07201229    
## 29 Mass Independant log(RMR)    0.0000000 0.13398490  0.18729777    
## 30 Mass Independant log(RMR)    1.2602254 0.14177631  0.20537453    
## 31 Mass Independant log(RMR)    1.6073708 0.15042976  0.21764843    
## 32 Mass Independant log(RMR)    1.6474074 0.12869909  0.17308051    
## 33 Mass Independant log(RMR)    2.7321064 0.13274864  0.17071098    
## 34 Mass Independant log(RMR)    2.8878776 0.13578794  0.19007496    
## 35 Mass Independant log(RMR)    3.2743628 0.14463523  0.20289714    
## 36 Mass Independant log(RMR)    4.0588194 0.14000948  0.18880081    
## 37 Mass Independant log(RMR)    4.0950899 0.09988873  0.21612058    
## 38 Mass Independant log(RMR)    4.5640317 0.14806374  0.20017371    
## 39 Mass Independant log(RMR)    4.8488012 0.13927880  0.19821267    
## 40 Mass Independant log(MMR)    0.0000000 0.36920443  0.37186427    
## 41 Mass Independant log(MMR)    0.5548350 0.37942598  0.38649376    
## 42 Mass Independant log(MMR)    0.6787176 0.36963415  0.36963415    
## 43 Mass Independant log(MMR)    0.7216353 0.37510675  0.37928176    
## 44 Mass Independant log(MMR)    0.8325546 0.37557342  0.37704611    
## 45 Mass Independant log(MMR)    1.2075921 0.38568092  0.39455042    
## 46 Mass Independant log(MMR)    1.3053415 0.37389091  0.37398454    
## 47 Mass Independant log(MMR)    1.4075558 0.38019832  0.38260583    
## 48 Mass Independant log(MMR)    2.6399397 0.37676558  0.37901010    
## 49 Mass Independant log(MMR)    3.2056969 0.38148663  0.38472604    
## 50 Mass Independant log(MMR)    3.2512482 0.37444184  0.37444184    
## 51 Mass Independant log(MMR)    3.3253542 0.38543858  0.39392723    
## 52 Mass Independant log(MMR)    3.3308469 0.38082566  0.38300719    
## 53 Mass Independant log(MMR)    3.5169143 0.38027035  0.38224935    
## 54 Mass Independant AAS         0.0000000 0.37930675  0.37930675    
## 55 Mass Independant AAS         1.5094369 0.38174466  0.38174466    
## 56 Mass Independant AAS         2.1236991 0.37949846  0.37949846    
## 57 Mass Independant AAS         2.6531066 0.38559940  0.38559940    
## 58 Mass Independant AAS         3.6519021 0.38196577  0.38196577    
## 59 Mass Independant AAS         4.3197068 0.37952426  0.37952426    
## 60 Mass Independant AAS         4.8437724 0.38574420  0.38574420    
## 61 Mass Independant CTM         0.0000000 0.30708516  0.30885939    
## 62 Mass Independant CTM         1.8460872 0.30784839  0.30858866    
##    Equation                                                                                                                                                                                                   
## 1  Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                             
## 2  Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                                          
## 3  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                
## 4  Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                             
## 5  Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                                                                                  
## 6  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                                                                                        
## 7  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                                                                     
## 8  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                                                                     
## 9  Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                                                                                  
## 10 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                             
## 11 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                
## 12 Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                                          
## 13 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                                                                                  
## 14 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                             
## 15 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                                                                     
## 16 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                                                                     
## 17 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                                                                                        
## 18 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                                                                                  
## 19 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc + Mother.Acc:Offspring.Acc                                                                          
## 20 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc + Mother.Acc:Offspring.Acc                                                                             
## 21 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc + Father.Acc:Offspring.Acc                                                                             
## 22 Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                                          
## 23 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                             
## 24 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                             
## 25 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                                                                                  
## 26 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                
## 27 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                                                                                  
## 28 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                                                                     
## 29 log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                                 
## 30 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                 
## 31 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                                                                      
## 32 Father.Acc + log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                    
## 33 Father.Acc + log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                                                                                            
## 34 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                    
## 35 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                                                         
## 36 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                                                                            
## 37 log(Mass) + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                                              
## 38 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc + Mother.Acc:Offspring.Acc                                                                 
## 39 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                                                         
## 40 log(Mass) + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                                              
## 41 log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                              
## 42 Father.Acc + log(Mass) + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                                 
## 43 log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                                 
## 44 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                 
## 45 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                 
## 46 Father.Acc + log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                    
## 47 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                    
## 48 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                                                                      
## 49 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                                                         
## 50 Father.Acc + log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                                                                                            
## 51 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                                                                      
## 52 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                                                                            
## 53 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                                                         
## 54 log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                              
## 55 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                 
## 56 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                                 
## 57 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                                                                      
## 58 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                                                    
## 59 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                                                                      
## 60 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                                                         
## 61 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc + Father.Acc:Offspring.Acc + Mother.Acc:Offspring.Acc + Father.Acc:Mother.Acc:Offspring.Acc            
## 62 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc + Father.Acc:Offspring.Acc + Mother.Acc:Offspring.Acc + Father.Acc:Mother.Acc:Offspring.Acc
write.csv(AIC.Selection, "Model_Outputs_AICc.csv", row.names = FALSE)

##Mass Specific Model Assessment ###Testing Model Assumptions for Resting Metabolic Rate Models: Model 1.

RMR.MS.m1 = lmer(log(RMR.MS) ~ Mother.Acc + Offspring.Acc + (1|Mother.ID) + (1|Father.ID),  
         REML = FALSE, data = RMR.MS.Data)

summary(RMR.MS.m1)$varcor 
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.08590 
##  Father.ID (Intercept) 0.15246 
##  Residual              0.41617
# Interesting. Father ID explains greater variance in data.

    Res.Alone = ggplot(data = data.frame("Res" = residuals(RMR.MS.m1, type = "deviance")), 
        aes(x = 1:nrow(RMR.MS.Data), y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
        ylab("Deviance Residuals")

    Res.by.Fit = ggplot(data = data.frame("Res" = residuals(RMR.MS.m1), "Fit" = predict(RMR.MS.m1)), 
        aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
        ylab("Deviance Residuals")

    Res.by.Resp = ggplot(data = data.frame("Res" = residuals(RMR.MS.m1), "Resp" = log(RMR.MS.Data$RMR.MS)), 
        aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
        xlab("log Mass Specific Metabolic Rate") + ylab("Deviance Residuals")
        
    QQNorm.1 = ggplot(data.frame("Res" = residuals(RMR.MS.m1), "Fit" = predict(RMR.MS.m1)),
        aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm.1, nrow = 2, top = 
        "Mass Specific Resting Metabolic Rate Residuals")  

    ## All looking well! Checking qqplot confidence intervals.

       qqPlot(residuals(RMR.MS.m1, type = "deviance"), ylab = "Deviance Residuals")

##  48 170 
##  43 145
    ## Good.

    r.squaredGLMM(RMR.MS.m1)
##             R2m       R2c
## [1,] 0.09689341 0.2325854
# Very low marginal R2, though not surprising. 

    # Parameter-specific heteroskedasticity
    {
    DamT.Res = ggplot(data = RMR.MS.Data, aes(x = as.factor(Mother.Acc), 
                  y = residuals(RMR.MS.m1, "deviance"))) + my.theme + theme_bw() + 
              geom_boxplot(fill = "mediumorchid", aes(middle = mean(residuals(RMR.MS.m1, "deviance")))) + 
            xlab("Dam Acclimation Temperature") + ylab("Deviance Residuals")

    OffT.Res = ggplot(data = RMR.MS.Data, aes(x = as.factor(Offspring.Acc), 
                 y = residuals(RMR.MS.m1, "deviance"))) + my.theme + theme_bw() + 
             geom_boxplot(fill = "mediumseagreen", aes(middle = mean(residuals(RMR.MS.m1, "deviance")))) + 
             xlab("Offspring Acclimation Temperature") + ylab("Deviance Residuals")
    }

    grid.arrange(DamT.Res, OffT.Res, nrow = 1, 
    top = "Fixed Effects by Residuals: Boxplot Bars Represent Means")  

## Residuals heterogenous.

###Model Result Plot.

dwdata = data.frame(term = c(rownames(summary(RMR.MS.m1)$coefficients)),
                    estimate = c(summary(RMR.MS.m1)$coefficients[1:3]),
                    std.error = c(summary(RMR.MS.m1)$coefficients[4:6]),
                    statistic = c(summary(RMR.MS.m1)$coefficients[7:9]))

p.values = c()
for (i in 1:nrow(dwdata)){
  p.values[i] = pt(abs(dwdata$statistic[i]), df = 8, lower.tail = FALSE)
} 
dwdata$p.values = round(p.values, 3)

dwdata$term = c("(Intercept)", "Dam Acclimation Temperature",
                "Offspring Acclimation Temperature")

dwdata$estimate = dwdata$estimate*10
dwdata$std.error= dwdata$std.error*10

dwdata %>%
  dwplot(style = "distribution", show_intercept = FALSE, alpha = 0.5) + theme_bw() + theme(legend.position = "none") + my.theme.dw + 
  xlab("Coefficient") + ylab("") + geom_vline(xintercept = 0, colour = "black", linetype = 2) +
  scale_fill_manual(values = "darkorchid") + scale_colour_manual(values = "black")  + 
  scale_x_continuous(limits = c(-5,5), breaks = c(-5.0, -2.5, 0, 2.5, 5.0), label = c("-0.50", "-0.25", "0", "0.25", "0.50")) +
  annotate("text", x = -4.2, y = 1.4, label = paste("p = ", dwdata$p.value[3])) + 
  annotate("text", x = 4.6, y = 2.2, label = paste("p = ", dwdata$p.value[2]))

dwdata$estimate = dwdata$estimate/10
dwdata$std.error= dwdata$std.error/10

###Testing Model Assumptions for Resting Metabolic Rate Models: Model 2.

RMR.MS.m2 = lmer(log(RMR.MS) ~ Offspring.Acc + (1|Mother.ID) +(1|Father.ID), REML = FALSE, 
        data = RMR.MS.Data)

summary(RMR.MS.m2)$varcor 
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.13616 
##  Father.ID (Intercept) 0.15553 
##  Residual              0.41631
# Slightly biased toward paternal ID.

    Res.Alone.m2 = ggplot(data = data.frame("Res" = residuals(RMR.MS.m2, type = "deviance")), 
        aes(x = 1:nrow(RMR.MS.Data), y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
        ylab("Deviance Residuals")

    Res.by.Fit.m2 = ggplot(data = data.frame("Res" = residuals(RMR.MS.m2), "Fit" = predict(RMR.MS.m2)), 
        aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
        ylab("Deviance Residuals")

    Res.by.Resp.m2 = ggplot(data = data.frame("Res" = residuals(RMR.MS.m2), "Resp" = log(RMR.MS.Data$RMR.MS)), 
        aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
        xlab("log Mass Specific Metabolic Rate") + ylab("Deviance Residuals")
        
    QQNorm.1.m2 = ggplot(data.frame("Res" = residuals(RMR.MS.m2), "Fit" = predict(RMR.MS.m2)),
        aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone.m2, Res.by.Fit.m2, Res.by.Resp.m2, QQNorm.1.m2, nrow = 2, top = 
        "Mass Specific Resting Metabolic Rate Residuals")  

    ## Again, looking well, outside of slight deviation at high quantiles. Checking qqplot confidence intervals.

       qqPlot(residuals(RMR.MS.m2, type = "deviance"), ylab = "Deviance Residuals")

##  48 170 
##  43 145
    # Still within confidence intervals.

    r.squaredGLMM(RMR.MS.m2)
##             R2m       R2c
## [1,] 0.05317331 0.2404377
    # Remarkably low marginal R2.. Arguably meaningless.

    ## Residual variance by factor.

    OffT.Res.m2 = ggplot(data = RMR.MS.Data, aes(x = as.factor(Offspring.Acc), 
                 y = residuals(RMR.MS.m1, "deviance"))) + my.theme + 
             theme_bw() + geom_boxplot(fill = "darkseagreen", aes(middle =
             mean(residuals(RMR.MS.m2, "deviance")))) + xlab("Offspring Acclimation Temperature") + 
             ylab("Deviance Residuals")

    print(OffT.Res.m2)

    # No distinct trend.

    # Calculating p-values.

    dwdata.m2 = data.frame(term = c(rownames(summary(RMR.MS.m2)$coefficients)),
                       estimate = c(summary(RMR.MS.m2)$coefficients[1:2]),
                       std.error = c(summary(RMR.MS.m2)$coefficients[3:4]),
                       statistic = c(summary(RMR.MS.m2)$coefficients[5:6]))

p.values.m2 = c()
for (i in 1:nrow(dwdata.m2)){
  p.values.m2[i] = pt(abs(dwdata.m2$statistic[i]), df = 8, lower.tail = FALSE)
} 
dwdata.m2$p.values = round(p.values.m2, 3)

dwdata.m2$term = c("(Intercept)", "Offspring Acclimation Temperature")

###Testing Model Assumptions for Resting Metabolic Rate Models: Model 3.

RMR.MS.m3 = lmer(log(RMR.MS) ~ Father.Acc + Mother.Acc + Offspring.Acc + 
        (1 | Mother.ID) + (1 | Father.ID), REML = FALSE, data = RMR.MS.Data)

summary(RMR.MS.m3)$varcor 
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.086156
##  Father.ID (Intercept) 0.131477
##  Residual              0.416194
# Similar variance explained to model iterations above. 

    Res.Alone.m3 = ggplot(data = data.frame("Res" = residuals(RMR.MS.m3, type = "deviance")), 
        aes(x = 1:nrow(RMR.MS.Data), y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
        ylab("Deviance Residuals")

    Res.by.Fit.m3 = ggplot(data = data.frame("Res" = residuals(RMR.MS.m3), "Fit" = predict(RMR.MS.m3)), 
        aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
        ylab("Deviance Residuals")

    Res.by.Resp.m3 = ggplot(data = data.frame("Res" = residuals(RMR.MS.m3), "Resp" = log(RMR.MS.Data$RMR.MS)), 
        aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
        xlab("log Mass Specific Metabolic Rate") + ylab("Deviance Residuals")
        
    QQNorm.1.m3 = ggplot(data.frame("Res" = residuals(RMR.MS.m3), "Fit" = predict(RMR.MS.m3)),
        aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone.m3, Res.by.Fit.m3, Res.by.Resp.m3, QQNorm.1.m3, nrow = 2, top = 
        "Mass Specific Resting Metabolic Rate Residuals")  

# Well behaved.

    r.squaredGLMM(RMR.MS.m3)
##            R2m       R2c
## [1,] 0.1194998 0.2294206
    # Interesting. Highest residual variation explained by fixed effects alone, of all model 
    # iterations.

## Residual variance by factor.

    {
    SireT.Res.m3 = ggplot(data = RMR.MS.Data, aes(x = as.factor(Father.Acc), 
                  y = residuals(RMR.MS.m3, "deviance"))) + my.theme + theme_bw() + 
              geom_boxplot(fill = "gold2", aes(middle = mean(residuals(RMR.MS.m3, "deviance")))) + 
            xlab("Dam Acclimation Temperature") + ylab("Deviance Residuals")

    DamT.Res.m3 = ggplot(data = RMR.MS.Data, aes(x = as.factor(Mother.Acc), 
                  y = residuals(RMR.MS.m3, "deviance"))) + my.theme + theme_bw() + 
              geom_boxplot(fill = "mediumorchid", aes(middle = mean(residuals(RMR.MS.m3, "deviance")))) + 
            xlab("Dam Acclimation Temperature") + ylab("Deviance Residuals")

    OffT.Res.m3 = ggplot(data = RMR.MS.Data, aes(x = as.factor(Offspring.Acc), 
                 y = residuals(RMR.MS.m3, "deviance"))) + my.theme + theme_bw() + 
             geom_boxplot(fill = "mediumseagreen", aes(middle = mean(residuals(RMR.MS.m3, "deviance")))) + 
             xlab("Offspring Acclimation Temperature") + ylab("Deviance Residuals")
    }

    grid.arrange(SireT.Res.m3, DamT.Res.m3, OffT.Res.m3, nrow = 1, 
    top = "Fixed Effects by Residuals: Boxplot Bars Represent Means")  

## Subtle distinctions in spread, but not of concern.

# Model result plots

dwdata.m3 = data.frame(term = c(rownames(summary(RMR.MS.m3)$coefficients)),
                       estimate = c(summary(RMR.MS.m3)$coefficients[1:4]),
                       std.error = c(summary(RMR.MS.m3)$coefficients[5:8]),
                       statistic = c(summary(RMR.MS.m3)$coefficients[9:12]))

p.values.m3 = c()
for (i in 1:nrow(dwdata.m3)){
  p.values.m3[i] = pt(abs(dwdata.m3$statistic[i]), df = 8, lower.tail = FALSE)
} 

dwdata.m3$p.values = round(p.values.m3, 3)

dwdata.m3$term = c("(Intercept)", "Sire Acclimation Temperature", "Dam Acclimation Temperature",
                   "Offspring Acclimation Temperature")

dwdata.m3$estimate = dwdata.m3$estimate*10
dwdata.m3$std.error= dwdata.m3$std.error*10

dwdata.m3 %>%
  dwplot(style = "distribution", show_intercept = FALSE) + theme_bw() + theme(legend.position = "none") + my.theme.dw + 
  xlab("Coefficient") + ylab("") + geom_vline(xintercept = 0, colour = "black", linetype = 2) +
  scale_fill_manual(values = "cornflowerblue") + scale_colour_manual(values = "black")  + 
  scale_x_continuous(limits = c(-6.0,6.0), breaks = c(-6.0, -3.0, 0, 3.0, 6.0), label = c("-0.60", "-0.30", "0", "0.30", "0.60")) +
  annotate("text", x = -4.0, y = 1.4, label = paste("p = ", dwdata.m3$p.value[4])) + 
  annotate("text", x = 4.5, y = 2.3, label = paste("p = ", dwdata.m3$p.value[3])) + 
  annotate("text", x = 4.3, y = 3.2, label = paste("p = ", dwdata.m3$p.value[2]))

dwdata.m3$estimate = dwdata.m3$estimate/10
dwdata.m3$std.error= dwdata.m3$std.error/10

Overlaying Resting Metabolic Rate Model Plots

dwdata = dwdata %>% mutate(model = "Model 1")
dwdata.m2 = dwdata.m2 %>% mutate(model = "Model 2")
dwdata.m3 = dwdata.m3 %>% mutate(model = "Model 3")

RMR.MS.Models = rbind(dwdata, dwdata.m2, dwdata.m3)

RMR.MS.Models$estimate = RMR.MS.Models$estimate*10
RMR.MS.Models$std.error = RMR.MS.Models$std.error*10 

RMR.MS.Models %>%
  dwplot(style = "distribution", show_intercept = FALSE) + theme_bw() + 
        theme(legend.justification = c(-0.1,0), 
            legend.position = c(0,0.01)) + my.theme.dw + xlab("Coefficient") + 
        ylab("") + geom_vline(xintercept = 0, 
            colour = "black", linetype = 2) + scale_fill_manual(values = c("mediumseagreen", 
        "mediumorchid", "cornflowerblue")) + scale_colour_manual(values = c("black", 
        "black", "black"))  + scale_x_continuous(limits = c(-6.0,6.0), 
            breaks = c(-6.0, -3.0, 0, 3.0, 6.0), 
        label = c("-0.60", "-0.30", "0", "0.30", "0.60"))

ggsave("/home/joshk/git_repositories/BrookTrout_TG/Figures/Mass Specific Resting Metabolic Rate Model Comparisons.jpeg", dpi = 800, limitsize = TRUE)

RMR.MS.Models$estimate = RMR.MS.Models$estimate/10
RMR.MS.Models$std.error = RMR.MS.Models$std.error/10 

RMR.Fit.m1 = data.frame("Response" = c(log(RMR.MS.Data$RMR.MS)), "Fit" = c(predict(RMR.MS.m1)))
RMR.Fit.m2 = data.frame("Response" = c(log(RMR.MS.Data$RMR.MS)), "Fit" = c(predict(RMR.MS.m2)))
RMR.Fit.m3 = data.frame("Response" = c(log(RMR.MS.Data$RMR.MS)), "Fit" = c(predict(RMR.MS.m3)))
RMR.Fit.m1 = RMR.Fit.m1 %>% mutate(model = "Model 1")
RMR.Fit.m2 = RMR.Fit.m2 %>% mutate(model = "Model 2")
RMR.Fit.m3 = RMR.Fit.m3 %>% mutate(model = "Model 3")

RMR.Fit.All = rbind(RMR.Fit.m1, RMR.Fit.m2, RMR.Fit.m3)
RMR.Fits = ggplot(data = RMR.Fit.All, aes(x = Fit, y = Response, fill = model, colour = model,
           linetype = model)) + theme_bw() + my.theme + geom_point(pch=21, size = 3, alpha = 0.5) + 
           geom_smooth(method = "lm") + xlab("Predicted log Mass Specific Resting Metabolic Rate") + 
           ylab("Measured log Mass Specific \n Resting Metabolic Rate") + scale_fill_manual(values =
           c("mediumseagreen", "mediumorchid", "cornflowerblue")) + scale_colour_manual(values = c("black", 
            "black", "black")) + theme(legend.title=element_blank())

print(RMR.Fits)

ggsave("/home/joshk/git_repositories/BrookTrout_TG/Figures/Mass Specific Resting Metabolic Rate Model Fit Comparisons.jpeg", dpi = 800, limitsize = TRUE)

##Mass Specific Model Assessment ###Testing Model Assumptions for Maximum Metabolic Rate Models: Model 1.

## Assessing parameter-specific heterogeneity.

    {
    PatT.MMR = ggplot(data = MMR.MS.Data, aes(x = as.factor(Father.Acc), 
        y = log(MMR.MS))) + theme_bw() + 
        geom_boxplot(fill = "mediumorchid", aes(middle = mean(log(MMR.MS)))) + 
        xlab("Sire Acclimation Temperature") + ylab("Log Mass Specific \n Maximum Metabolic Rate \n (mg O2/kg/h)")

    DamT.MMR = ggplot(data = MMR.MS.Data, aes(x = as.factor(Mother.Acc), 
        y = log(MMR.MS))) + theme_bw() + 
        geom_boxplot(fill = "mediumseagreen", aes(middle = mean(log(MMR.MS)))) + 
        xlab("Dam Acclimation Temperature") + ylab("Log Mass Specific \n Maximum Metabolic Rate \n (mg O2/kg/h)")

    OffT.MMR = ggplot(data = MMR.MS.Data, aes(x = as.factor(Offspring.Acc), 
        y = log(MMR.MS))) + theme_bw() + 
        geom_boxplot(fill = "royalblue", aes(middle = mean(log(MMR.MS)))) + 
        xlab("Offspring Acclimation Temperature") + 
    ylab("Log Mass Specific \n Maximum Metabolic Rate \n (mg O2/kg/h)")

    Pat.by.MMR = ggplot(data = MMR.MS.Data, aes(x = as.factor(Father.Acc), 
        y = log(MMR.MS), fill = Offspring.Acc)) + theme_bw() + 
        theme(legend.position = "top") + geom_boxplot() + 
        scale_fill_manual(values = c("mediumorchid", "royalblue")) + 
        xlab("Sire Acclimation Temperature") + 
       ylab("Log Mass Specific \n Maximum Metabolic Rate \n (mg O2/kg/h)")

    Pat.by.MMR = ggplot(data = MMR.MS.Data, aes(x = as.factor(Father.Acc), 
        y = log(MMR.MS), fill = Mother.Acc)) + theme_bw() + 
        theme(legend.position = "top") + geom_boxplot() + scale_fill_manual(values = c("mediumseagreen", 
        "royalblue")) + xlab("Sire Acclimation Temperature") + 
      ylab("Log Mass Specific \n Maximum Metabolic Rate \n (mg O2/kg/h)")

    Off.by.MMR = ggplot(data = MMR.MS.Data, aes(x = as.factor(Offspring.Acc), 
        y = log(MMR.MS), fill = Mother.Acc)) + theme_bw() + 
        theme(legend.position = "top") + geom_boxplot() + scale_fill_manual(values = c("mediumorchid", 
        "mediumseagreen")) + xlab("Offspring Acclimation Temperature") + 
       ylab("Log Mass Specific \n Maximum Metabolic Rate \n (mg O2/kg/h)")

    }

    grid.arrange(PatT.MMR, DamT.MMR, OffT.MMR, Pat.by.MMR, Pat.by.MMR,
    Off.by.MMR, nrow = 2, 
    top = "Fixed Effects by Response: Boxplot Bars Represent Means") 

    require('plyr')

    {
    PatT.MMR.Err = ggplot(ddply(MMR.MS.Data, ~Father.Acc, summarise, mean = mean(log(MMR.MS)),
        High.Conf = mean(log(MMR.MS)) + 2*sd(log(MMR.MS)), Low.Conf = mean(log(MMR.MS)) - 
        2*sd(log(MMR.MS)))) + theme_bw() + geom_errorbar(mapping = aes(x = as.factor(Father.Acc), 
        ymin = c(Low.Conf), ymax = c(High.Conf)), width = 0.2, size = 1, color = "cornflowerblue") + 
        xlab("Sire Acclimation Temperature") + ylab("Log Mass Specific \n Maximum Metabolic Rate \n (mg O2/kg/h)")
        
    DamT.MMR.Err = ggplot(ddply(MMR.MS.Data, ~Mother.Acc, summarise, mean = mean(log(MMR.MS)),
        High.Conf = mean(log(MMR.MS)) + 2*sd(log(MMR.MS)), Low.Conf = mean(log(MMR.MS)) - 
        2*sd(log(MMR.MS)))) + theme_bw() + geom_errorbar(mapping = aes(x = as.factor(Mother.Acc), 
        ymin = c(Low.Conf), ymax = c(High.Conf)), width = 0.2, size = 1, color = "mediumorchid") + 
        xlab("Dam Acclimation Temperature") + ylab("Log Mass Specific \n Maximum Metabolic Rate \n (mg O2/kg/h)")

    OffT.MMR.Err = ggplot(ddply(MMR.MS.Data, ~Offspring.Acc, summarise, mean = mean(log(MMR.MS)),
        High.Conf = mean(log(MMR.MS)) + 2*sd(log(MMR.MS)), Low.Conf = mean(log(MMR.MS)) - 
        2*sd(log(MMR.MS)))) + theme_bw() + geom_errorbar(mapping = aes(x = as.factor(Offspring.Acc), 
        ymin = c(Low.Conf), ymax = c(High.Conf)), width = 0.2, size = 1, color = "mediumseagreen") + 
        xlab("Offspring Acclimation Temperature") + 
      ylab("Log Mass Specific \n Maximum Metabolic Rate \n (mg O2/kg/h)")

    Pat.by.Off.MMR.Err = ggplot(ddply(MMR.MS.Data, ~Father.Acc*Offspring.Acc, summarise, 
                                      mean = mean(log(MMR.MS)),
        High.Conf = mean(log(MMR.MS)) + 2*sd(log(MMR.MS)), Low.Conf = mean(log(MMR.MS)) - 
        2*sd(log(MMR.MS))), aes(colour = Offspring.Acc)) + theme_bw() + 
      geom_errorbar(mapping = aes(x = as.factor(Father.Acc), 
        ymin = c(Low.Conf), ymax = c(High.Conf)), width = 0.2, size = 1) + 
        scale_colour_manual(values = c("cornflowerblue", "gold2")) + theme(legend.position = "top") + 
        xlab("Sire Acclimation Temperature") + ylab("Log Mass Specific \n Maximum Metabolic Rate \n (mg O2/kg/h)")

    Pat.by.Mat.MMR.Err = ggplot(ddply(MMR.MS.Data, ~Father.Acc*Mother.Acc, summarise, mean = mean(log(MMR.MS)),
        High.Conf = mean(log(MMR.MS)) + 2*sd(log(MMR.MS)), Low.Conf = mean(log(MMR.MS)) - 
        2*sd(log(MMR.MS))), aes(colour = Mother.Acc)) + theme_bw() + 
      geom_errorbar(mapping = aes(x = as.factor(Father.Acc), 
        ymin = c(Low.Conf), ymax = c(High.Conf)), width = 0.2, size = 1) + 
        scale_colour_manual(values = c("cornflowerblue", "violetred2")) + theme(legend.position = "top") +
        xlab("Sire Acclimation Temperature") + ylab("Log Mass Specific \n Maximum Metabolic Rate \n (mg O2/kg/h)")

    Mat.by.Off.MMR.Err = ggplot(ddply(MMR.MS.Data, ~Offspring.Acc*Mother.Acc, summarise, 
                                      mean = mean(log(MMR.MS)),
        High.Conf = mean(log(MMR.MS)) + 2*sd(log(MMR.MS)), Low.Conf = mean(log(MMR.MS)) - 
        2*sd(log(MMR.MS))), aes(colour = Offspring.Acc)) + theme_bw() + 
      geom_errorbar(mapping = aes(x = as.factor(Mother.Acc), 
        ymin = c(Low.Conf), ymax = c(High.Conf)), width = 0.2, size = 1) + 
        scale_colour_manual(values = c("gold2", "violetred2")) + theme(legend.position = "top") +
        xlab("Dam Acclimation Temperature") + ylab("Log Mass Specific \n Maximum Metabolic Rate \n (mg O2/kg/h)")
    }

    grid.arrange(PatT.MMR.Err, DamT.MMR.Err, OffT.MMR.Err, Pat.by.Off.MMR.Err, 
            Pat.by.Mat.MMR.Err, Mat.by.Off.MMR.Err, nrow = 2, 
    top = "Fixed Effects by Response: Error Bars Represent 95% Confidence Intervals") 

    # Little reason for concern.

## Model 1

MMR.MS.m1 = lmer(log(MMR.MS) ~ Father.Acc + Offspring.Acc + 
            (1|Mother.ID) + (1|Father.ID), REML = FALSE, data = MMR.MS.Data)

summary(MMR.MS.m1)$varcor 
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.00000 
##  Father.ID (Intercept) 0.03177 
##  Residual              0.20213
# Maternal ID explains zero variance. Retaining for comparability with downstream model iterations.

    Res.Alone.m1 = ggplot(data = data.frame("Res" = residuals(MMR.MS.m1, type = "deviance")), 
        aes(x = 1:nrow(MMR.MS.Data), y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
        ylab("Deviance Residuals")

    Res.by.Fit.m1 = ggplot(data = data.frame("Res" = residuals(MMR.MS.m1), "Fit" = predict(MMR.MS.m1)), 
        aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
        ylab("Deviance Residuals")

    Res.by.Resp.m1 = ggplot(data = data.frame("Res" = residuals(MMR.MS.m1), "Resp" = log(MMR.MS.Data$MMR.MS)), 
        aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
        xlab("log Mass Specific Metabolic Rate") + ylab("Deviance Residuals")
        
    QQNorm.1.m1 = ggplot(data.frame("Res" = residuals(MMR.MS.m1), "Fit" = predict(MMR.MS.m1)),
        aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone.m1, Res.by.Fit.m1, Res.by.Resp.m1, QQNorm.1.m1, nrow = 2, top = 
        "Mass Specific Maximum Metabolic Rate Residuals")

# Residuals vs fitted okay, but residuals very non-normal.

    ## Identifying individuals of concern.
    
    qqPlot(residuals(MMR.MS.m1, type = "deviance"), ylab = "Deviance Residuals")

## 170  86 
## 167  84
    # Majority of right tail exceeding confidence intervals.

    {
    PatT.Res.m1 = ggplot(data = MMR.MS.Data, aes(x = as.factor(Father.Acc), 
        y = residuals(MMR.MS.m1, "deviance"))) + theme_bw() + 
        geom_boxplot(fill = "mediumorchid", aes(middle = mean(residuals(MMR.MS.m1, "deviance")))) + 
        xlab("Sire Acclimation Temperature") + ylab("Deviance Residuals")

    OffT.Res.m1 = ggplot(data = MMR.MS.Data, aes(x = as.factor(Offspring.Acc), 
        y = residuals(MMR.MS.m1, "deviance"))) + theme_bw() + 
        geom_boxplot(fill = "royalblue", aes(middle = mean(residuals(MMR.MS.m1, "deviance")))) + 
        xlab("Offspring Acclimation Temperature") + ylab("Deviance Residuals")

    Pat.by.Off.T.Res.m1 = ggplot(data = MMR.MS.Data, aes(x = as.factor(Father.Acc), 
        y = residuals(MMR.MS.m1, type = "deviance"), fill = Offspring.Acc)) + theme_bw() + 
        theme(legend.position = "top") + geom_boxplot() + 
        scale_fill_manual(values = c("mediumorchid", "royalblue")) + 
        xlab("Sire Acclimation Temperature") + ylab("Deviance Residuals")

    }
    grid.arrange(PatT.Res.m1, OffT.Res.m1, Pat.by.Off.T.Res.m1, nrow = 2, 
    top = "Fixed Effects by Residuals: Boxplot Bars Represent Means") 

    # No concern with variance structure.

    # Assessing outlier patterns.
    
    Out = MMR.MS.Data[c(which(residuals(MMR.MS.m1, type = "deviance") >= 0.22)),]
    In = MMR.MS.Data[-c(which(residuals(MMR.MS.m1, type = "deviance") >= 0.22)),]

    Out = Out %>% mutate(Group = "Outlier")
    In = In %>% mutate(Group = "Inlier")
    MMR.MS.Data.mod = rbind(Out, In)

        Out.Mass.MMR = ggplot(data = MMR.MS.Data.mod, aes(x = as.factor(Group), 
        y = Mass)) + theme_bw() + 
        geom_boxplot(fill = "royalblue") + 
        xlab("Outlier Group") + ylab("Mass (g)")

        Out.MMR = ggplot(data = MMR.MS.Data.mod, aes(x = as.factor(Group), 
        y = MMR.MS)) + theme_bw() + 
        geom_boxplot(fill = "royalblue") + 
        xlab("Outlier Group") + ylab("Maximum Mass Specific \n Metabolic Rate")

    grid.arrange(Out.Mass.MMR, Out.MMR)

    Sweep = as.data.frame(ddply(MMR.MS.Data.mod, ~Group, summarise, 
        Off.Ratio = sum(Offspring.Acc == 1)/sum(Offspring.Acc == 2),
        Pat.Ratio = sum(Father.Acc == 1)/sum(Father.Acc == 2)))

    Global = data.frame("Group" = "Global", "Off.Ratio" = length(which(MMR.MS.Data.mod$Offspring.Acc ==1))/
            length(which(MMR.MS.Data.mod$Offspring.Acc ==2)), 
            "Pat.Ratio" = length(which(MMR.MS.Data.mod$Father.Acc ==1))/
            length(which(MMR.MS.Data.mod$Father.Acc ==2)))

    Sweep = rbind(Sweep, Global)
    
    # Disproportionate number of individuals from cold-reared fathers in outliers.
        
    # Perhaps test a constant variance by sire acclimation temperature.

    MMR.MS.m1 = lme(log(MMR.MS) ~ Father.Acc + Offspring.Acc, random =  
            list(~1|Mother.ID, ~1|Father.ID), weights = varIdent(form = ~1|Father.Acc), 
        data = MMR.MS.Data)

    qqPlot(residuals(MMR.MS.m1, type = "normalized"), ylab = "Normalized Residuals")

## 3/4 1/3 
##  84 167
    # Very little, if any, correction.

    MMR.MS.m1 = lmer(log(MMR.MS) ~ Father.Acc + Offspring.Acc + 
            (1|Mother.ID) + (1|Father.ID), REML = FALSE, data = MMR.MS.Data)

    Resid.hist = ggplot(data.frame(Res = c(residuals(MMR.MS.m1, type = "deviance")), Model = "Model1"), 
    aes(Res, fill = Model, colour = Model)) + geom_histogram(alpha = 0.5) + theme_bw() + 
    my.theme + xlab("Deviance Residuals") + ylab("Count") + scale_fill_manual(values = 
    c("cornflowerblue")) + scale_colour_manual(values = c("black")) + theme(legend.position = "none") + 
    geom_density(alpha = 0.5, aes(y=0.3/8 * ..count..))

    Resp.hist = ggplot(data.frame(Resp = c(log(MMR.MS.Data$MMR.MS)), Model = "Model1"), 
    aes(Resp, fill = Model, colour = Model)) + geom_histogram(alpha = 0.5) + theme_bw() + 
    my.theme + xlab("log Mass Specific Maximum \n Metabolic Rate") + ylab("Count") + scale_fill_manual(values = 
    c("gold2")) + scale_colour_manual(values = c("black")) + theme(legend.position = "none") + 
    geom_density(alpha = 0.5, aes(y=0.25/7.5 * ..count..))

    grid.arrange(Resid.hist, Resp.hist, nrow = 1, top = 
        "Mass Specific Maximum Metabolic Rate Histograms")

    ## Strangely bimodal, but *no more bimodal than the raw response*. Perhaps improper to use these data at all?

# Calculating p-values

MMR.dwdata = data.frame(term = c(rownames(summary(MMR.MS.m1)$coefficients)),
                         estimate = c(summary(MMR.MS.m1)$coefficients[1:3]),
                         std.error = c(summary(MMR.MS.m1)$coefficients[4:6]),
                         statistic = c(summary(MMR.MS.m1)$coefficients[7:9]))

p.values = c()
for (i in 1:nrow(MMR.dwdata)){
  p.values[i] = pt(abs(MMR.dwdata$statistic[i]), df = 8, lower.tail = FALSE)
} 
MMR.dwdata$p.values = round(p.values, 3)

MMR.dwdata$term = c("(Intercept)", "Sire Acclimation Temperature", 
    "Offspring Acclimation Temperature")

print(MMR.dwdata) 
##                                term   estimate  std.error  statistic p.values
## 1                       (Intercept)  5.9923085 0.03258909 183.874659    0.000
## 2      Sire Acclimation Temperature  0.1027712 0.04174533   2.461860    0.020
## 3 Offspring Acclimation Temperature -0.1181378 0.02724302  -4.336442    0.001

###Testing Model Assumptions for Maximum Metabolic Rate Models: Model 2.

MMR.MS.m2 = lmer(log(MMR.MS) ~ Father.Acc + Mother.Acc + Offspring.Acc + 
        (1 | Mother.ID) + (1 | Father.ID), REML = FALSE, data = MMR.MS.Data)

summary(MMR.MS.m2)$varcor 
##  Groups    Name        Std.Dev.  
##  Mother.ID (Intercept) 3.7095e-06
##  Father.ID (Intercept) 3.2569e-02
##  Residual              2.0162e-01
# Both mother and dather ID explain minimal, but recognizable variance.

    Res.Alone.m2 = ggplot(data = data.frame("Res" = residuals(MMR.MS.m2, type = "deviance")), 
        aes(x = 1:nrow(MMR.MS.Data), y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
        ylab("Deviance Residuals")

    Res.by.Fit.m2 = ggplot(data = data.frame("Res" = residuals(MMR.MS.m2), "Fit" = predict(MMR.MS.m2)), 
        aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
        ylab("Deviance Residuals")

    Res.by.Resp.m2 = ggplot(data = data.frame("Res" = residuals(MMR.MS.m2), "Resp" = log(MMR.MS.Data$MMR.MS)), 
        aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
        xlab("log Mass Specific Metabolic Rate") + ylab("Deviance Residuals")
        
    QQNorm.1.m2 = ggplot(data.frame("Res" = residuals(MMR.MS.m2), "Fit" = predict(MMR.MS.m2)),
        aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone.m2, Res.by.Fit.m2, Res.by.Resp.m2, QQNorm.1.m2, nrow = 2, top = 
        "Mass Specific Resting Metabolic Rate Residuals")

## Similar trends to model 1.

    qqPlot(residuals(MMR.MS.m2, type = "deviance"), ylab = "Deviance Residuals")

## 170  86 
## 167  84
    Resid.hist = ggplot(data.frame(Res = c(residuals(MMR.MS.m2, type = "deviance")), Model = "Model1"), 
    aes(Res, fill = Model, colour = Model)) + geom_histogram(alpha = 0.5) + theme_bw() + 
    my.theme + xlab("Deviance Residuals") + ylab("Count") + scale_fill_manual(values = 
    c("cornflowerblue")) + scale_colour_manual(values = c("black")) + theme(legend.position = "none") + 
    geom_density(alpha = 0.5, aes(y=0.3/8 * ..count..))

    print(Resid.hist)

    # Bimodal as before.

 #Parameter-specific heterogeneity

    {
    PatT.Res.m2 = ggplot(data = MMR.MS.Data, aes(x = as.factor(Father.Acc), 
        y = residuals(MMR.MS.m2, "deviance"))) + theme_bw() + 
        geom_boxplot(fill = "mediumorchid", aes(middle = mean(residuals(MMR.MS.m2, "deviance")))) + 
        xlab("Sire Acclimation Temperature") + ylab("Deviance Residuals")

    OffT.Res.m2 = ggplot(data = MMR.MS.Data, aes(x = as.factor(Offspring.Acc), 
        y = residuals(MMR.MS.m2, "deviance"))) + theme_bw() + 
        geom_boxplot(fill = "royalblue", aes(middle = mean(residuals(MMR.MS.m2, "deviance")))) + 
        xlab("Offspring Acclimation Temperature") + ylab("Deviance Residuals")

    DamT.Res.m2 = ggplot(data = MMR.MS.Data, aes(x = as.factor(Mother.Acc), 
        y = residuals(MMR.MS.m2, type = "deviance"))) + theme_bw() + 
        geom_boxplot(fill = "royalblue", aes(middle = mean(residuals(MMR.MS.m2, "deviance")))) + 
        scale_fill_manual(values = c("mediumorchid", "royalblue")) + 
        xlab("Sire Acclimation Temperature") + ylab("Deviance Residuals")

    }
    grid.arrange(PatT.Res.m2, OffT.Res.m2, DamT.Res.m2, nrow = 2, 
    top = "Fixed Effects by Residuals: Boxplot Bars Represent Means") 

    # As above, nothing of concern.

r.squaredGLMM(MMR.MS.m2)
##            R2m       R2c
## [1,] 0.1248206 0.1470761
# Expectedly low.

# Model result plots

MMR.dwdata2 = data.frame(term = c(rownames(summary(MMR.MS.m2)$coefficients)),
                         estimate = c(summary(MMR.MS.m2)$coefficients[1:4]),
                         std.error = c(summary(MMR.MS.m2)$coefficients[5:8]),
                         statistic = c(summary(MMR.MS.m2)$coefficients[9:12]))

p.values = c()
for (i in 1:nrow(MMR.dwdata2)){
  p.values[i] = pt(abs(MMR.dwdata2$statistic[i]), df = 8, lower.tail = FALSE)
} 

MMR.dwdata2$p.values = round(p.values, 3)

MMR.dwdata2$term = c("(Intercept)", "Sire Acclimation Temperature", 
        "Dam Acclimation Temperature", "Offspring Acclimation Temperature")

print(MMR.dwdata2) 
##                                term    estimate  std.error  statistic p.values
## 1                       (Intercept)  5.97655289 0.03652341 163.636215    0.000
## 2      Sire Acclimation Temperature  0.10402494 0.04233018   2.457465    0.020
## 3       Dam Acclimation Temperature  0.02715208 0.02712690   1.000928    0.173
## 4 Offspring Acclimation Temperature -0.11767031 0.02718229  -4.328933    0.001
MMR.dwdata2$estimate = MMR.dwdata2$estimate*20
MMR.dwdata2$std.error= MMR.dwdata2$std.error*20

MMR.dwdata2 %>%
  dwplot(style = "distribution", show_intercept = FALSE) + theme_bw() + theme(legend.position = "none") + my.theme.dw + 
  xlab("Coefficient") + ylab("") + geom_vline(xintercept = 0, colour = "black", linetype = 2) +
  scale_fill_manual(values = "cornflowerblue") + scale_colour_manual(values = "black")  + 
  scale_x_continuous(limits = c(-5.0,5.0), breaks = c(-5.0, -2.5, 0, 2.5, 5.0), label = c("-0.25", "-0.13", "0", "0.13", "0.25")) +
  annotate("text", x = -4.0, y = 1.4, label = paste("p = ", MMR.dwdata2$p.value[4])) + 
  annotate("text", x = 2.5, y = 2.3, label = paste("p = ", MMR.dwdata2$p.value[3])) + 
  annotate("text", x = 4.5, y = 3.2, label = paste("p = ", MMR.dwdata2$p.value[2]))

MMR.dwdata2$estimate = MMR.dwdata2$estimate/20
MMR.dwdata2$std.error= MMR.dwdata2$std.error/20

###Testing Model Assumptions for Maximum Metabolic Rate Models: Model 3.

MMR.MS.m3 = lmer(log(MMR.MS) ~ Offspring.Acc +  
            (1|Mother.ID) + (1|Father.ID), REML = FALSE, data = MMR.MS.Data)

summary(MMR.MS.m3)$varcor 
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.000000
##  Father.ID (Intercept) 0.060574
##  Residual              0.202140
# Similar to model 1, maternal ID explains zero variance.

    Res.Alone.m3 = ggplot(data = data.frame("Res" = residuals(MMR.MS.m3, type = "deviance")), 
        aes(x = 1:nrow(MMR.MS.Data), y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
        ylab("Deviance Residuals")

    Res.by.Fit.m3 = ggplot(data = data.frame("Res" = residuals(MMR.MS.m3), "Fit" = predict(MMR.MS.m3)), 
        aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
        ylab("Deviance Residuals")

    Res.by.Resp.m3 = ggplot(data = data.frame("Res" = residuals(MMR.MS.m3), "Resp" = log(MMR.MS.Data$MMR.MS)), 
        aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
        xlab("log Mass Specific Metabolic Rate") + ylab("Deviance Residuals")
        
    QQNorm.1.m3 = ggplot(data.frame("Res" = residuals(MMR.MS.m3), "Fit" = predict(MMR.MS.m3)),
        aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone.m3, Res.by.Fit.m3, Res.by.Resp.m3, QQNorm.1.m3, nrow = 2, top = 
        "Mass Specific Resting Metabolic Rate Residuals")

## Right residual tail persists.

    qqPlot(residuals(MMR.MS.m3, type = "deviance"), ylab = "Deviance Residuals")

## 170  86 
## 167  84
r.squaredGLMM(MMR.MS.m3)
##             R2m      R2c
## [1,] 0.07304761 0.149428
# Histogram of residuals?

ggplot(data.frame(Res = c(residuals(MMR.MS.m3, type = "deviance")), Model = "Model1"), 
    aes(Res, fill = Model, colour = Model)) + geom_histogram(alpha = 0.5) + theme_bw() + 
    my.theme + xlab("Deviance Residuals") + ylab("Count") + scale_fill_manual(values = 
    c("cornflowerblue")) + scale_colour_manual(values = c("black")) + theme(legend.position = "none") + 
    geom_density(alpha = 0.5, aes(y=0.3/8 * ..count..))

# Further subtle bimodality.

MMR.Off.Hist = ggplot(data.frame(Res = c(residuals(MMR.MS.m3, type = "deviance")), 
                                 Off.Acc = MMR.MS.Data$Offspring.Acc),
    aes(Res, fill = Off.Acc, colour = Off.Acc)) + geom_histogram(alpha = 0.5, bins = 20) + theme_bw() + 
    my.theme + xlab("Pearson Residuals") + ylab("Count") + scale_fill_manual(values = 
    c("mediumorchid", "royalblue")) + scale_colour_manual(values = c("black", "black")) + 
    geom_vline(xintercept = mean(residuals(MMR.MS.m3, type = "deviance")), linetype="dashed", size = 1) + 
    theme(legend.position = "top")

MMR.Pat.Hist = ggplot(data.frame(Res = c(residuals(MMR.MS.m3, type = "deviance")), 
                                 Pat.Acc = MMR.MS.Data$Father.Acc),
    aes(Res, fill = Pat.Acc, colour = Pat.Acc)) + geom_histogram(alpha = 0.5, bins = 20) + theme_bw() + 
    my.theme + xlab("Pearson Residuals") + ylab("Count") + scale_fill_manual(values = 
    c("mediumorchid", "royalblue")) + scale_colour_manual(values = c("black", "black")) + 
    geom_vline(xintercept = mean(residuals(MMR.MS.m3, type = "deviance")), linetype="dashed", size = 1) + 
    theme(legend.position = "top")

MMR.PG.Hist = ggplot(data.frame(Res = c(residuals(MMR.MS.m3, type = "deviance")), 
                                PG = as.factor(MMR.MS.Data$Parent.group)),
    aes(Res, fill = PG, colour = PG)) + geom_histogram(alpha = 0.5, bins = 20) + theme_bw() + 
    my.theme + xlab("Pearson Residuals") + ylab("Count") + scale_fill_manual(values = 
    c("mediumorchid", "royalblue", "mediumseagreen", "gold2")) + 
    scale_colour_manual(values = c("black", "black", "black", "black")) + 
    geom_vline(xintercept = mean(residuals(MMR.MS.m3, type = "deviance")), linetype="dashed", size = 1) + 
    theme(legend.position = "top")

    grid.arrange(MMR.Off.Hist, MMR.Pat.Hist, MMR.PG.Hist, nrow = 1, top = 
        "Mass Specific Maximum Metabolic Rate Residuals Histograms")

## Oddly bimodal accross all groups, though seemingly worse in cold father treatment.

## Exploring parameter-specific heteroskedasticity

    MMR.All.Res = ggplot(data = MMR.MS.Data, aes(x = as.factor(Offspring.Acc), 
        y = residuals(MMR.MS.m3, "deviance"), fill = as.factor(Offspring.Acc))) + 
        my.theme + theme_bw() + geom_boxplot(lwd=1) + scale_fill_manual(values=c("ivory", 
        "royalblue3")) + 
        xlab("Sire Acclimation Temperature") + ylab("Deviance Residuals") + 
        theme(legend.position="top")

# Residuals appear evenly dispersed between groups 
    
# Model result plots

MMR.dwdata3 = data.frame(term = c(rownames(summary(MMR.MS.m3)$coefficients)),
                         estimate = c(summary(MMR.MS.m3)$coefficients[1:2]),
                         std.error = c(summary(MMR.MS.m3)$coefficients[3:4]),
                         statistic = c(summary(MMR.MS.m3)$coefficients[5:6]))

p.values = c()
for (i in 1:nrow(MMR.dwdata3)){
  p.values[i] = pt(abs(MMR.dwdata3$statistic[i]), df = 8, lower.tail = FALSE)
} 

MMR.dwdata3$p.values = round(p.values, 3)

MMR.dwdata3$term = c("(Intercept)", "Offspring Acclimation Temperature")

MMR.dwdata3$estimate = MMR.dwdata3$estimate*20
MMR.dwdata3$std.error= MMR.dwdata3$std.error*20

MMR.dwdata3 %>%
  dwplot(style = "distribution", show_intercept = FALSE) + theme_bw() + theme(legend.position = "none") +
  my.theme.dw + 
  xlab("Coefficient") + ylab("") + geom_vline(xintercept = 0, colour = "black", linetype = 2) +
  scale_fill_manual(values = "mediumseagreen") + scale_colour_manual(values = "black")  + 
  scale_x_continuous(limits = c(-6.0,6.0), breaks = c(-6.0, -3.0, 0, 3.0, 6.0), 
                     label = c("-0.30", "-0.15","0", "0.15", "0.30")) +
  annotate("text", x = -5.0, y = 1.5, label = paste("p =\n", MMR.dwdata3$p.value[2]))

MMR.dwdata3$estimate = MMR.dwdata3$estimate/20
MMR.dwdata3$std.error= MMR.dwdata3$std.error/20

Overlaying MMR model plots

MMR.dwdata = MMR.dwdata %>% mutate(model = "Model 1")
MMR.dwdata2 = MMR.dwdata2 %>% mutate(model = "Model 2")
MMR.dwdata3 = MMR.dwdata3 %>% mutate(model = "Model 3")

MMR.MS.Models = rbind(MMR.dwdata, MMR.dwdata2, MMR.dwdata3)

MMR.MS.Models$estimate = MMR.MS.Models$estimate*20
MMR.MS.Models$std.error = MMR.MS.Models$std.error*20 

MMR.MS.Models %>%
  dwplot(style = "distribution", show_intercept = FALSE) + theme_bw() + 
  theme(legend.justification = c(-6.0,-0.1), 
  legend.position = c(0,0.01)) + my.theme.dw + xlab("Coefficient") + ylab("") + geom_vline(xintercept = 0, 
  colour = "black", linetype = 2) + scale_fill_manual(values = 
  c("mediumseagreen", "mediumorchid", "cornflowerblue")) + 
  scale_colour_manual(values = c("black", "black", "black"))  + scale_x_continuous(limits = c(-6.0,6.0), 
  breaks = c(-6.0, -3.0, 0, 3.0, 6.0), label = c("-0.30", "-0.15", "0", "0.15", "0.30"))

ggsave("/home/joshk/git_repositories/BrookTrout_TG/Figures/Mass Specific Maximum Metabolic Rate Model Comparisons.jpeg", dpi = 800, limitsize = TRUE)

MMR.MS.Models$estimate = MMR.MS.Models$estimate/20
MMR.MS.Models$std.error = MMR.MS.Models$std.error/20 

MMR.Fit.m1 = data.frame("Response" = c(log(MMR.MS.Data$MMR.MS)), "Fit" = c(predict(MMR.MS.m1)))
MMR.Fit.m2 = data.frame("Response" = c(log(MMR.MS.Data$MMR.MS)), "Fit" = c(predict(MMR.MS.m2)))
MMR.Fit.m3 = data.frame("Response" = c(log(MMR.MS.Data$MMR.MS)), "Fit" = c(predict(MMR.MS.m3)))
MMR.Fit.m1 = MMR.Fit.m1 %>% mutate(model = "Model 1")
MMR.Fit.m2 = MMR.Fit.m2 %>% mutate(model = "Model 2")
MMR.Fit.m3 = MMR.Fit.m3 %>% mutate(model = "Model 3")

MMR.Fit.All = rbind(MMR.Fit.m1, MMR.Fit.m2, MMR.Fit.m3)
MMR.Fits = ggplot(data = MMR.Fit.All, aes(x = Fit, y = Response, fill = model, colour = model,
           linetype = model)) + theme_bw() + my.theme + geom_point(pch=21, size = 3, alpha = 0.5) + 
           geom_smooth(method = "lm") + xlab("Predicted log Mass Specific Maximum Metabolic Rate") + 
           ylab("Measured log Mass Specific \n Maximum Metabolic Rate") + 
           scale_fill_manual(values = c("mediumseagreen", 
           "mediumorchid", "cornflowerblue")) + scale_colour_manual(values = c("black", 
         "black", "black")) + theme(legend.title=element_blank())

print(MMR.Fits)

ggsave("/home/joshk/git_repositories/BrookTrout_TG/Figures/log Mass Specific Maximum Metabolic Rate Fit Comparisons - AICc.jpeg", dpi = 800, limitsize = TRUE)

##Mass Specific Model Assessment ###Testing Model Assumptions for Absolute Aerobic Scope Models: Model 1.

AAS.MS.m1 = lmer(log(AAS.MS) ~ Offspring.Acc + (1|Mother.ID) + 
                   (1|Father.ID), REML = FALSE, data = AAS.MS.Data)

summary(AAS.MS.m1)$varcor 
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.00000 
##  Father.ID (Intercept) 0.00000 
##  Residual              0.29537
# Interestign. Zero variance explained by either random intercept.

    Res.Alone.m1 = ggplot(data = data.frame("Res" = residuals(AAS.MS.m1, type = "deviance")), 
        aes(x = 1:nrow(AAS.MS.Data), y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
        ylab("Deviance Residuals")

    Res.by.Fit.m1 = ggplot(data = data.frame("Res" = residuals(AAS.MS.m1), "Fit" = predict(AAS.MS.m1)), 
        aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
        ylab("Deviance Residuals")

    Res.by.Resp.m1 = ggplot(data = data.frame("Res" = residuals(AAS.MS.m1), "Resp" = log(AAS.MS.Data$AAS.MS)), 
        aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
        xlab("log Mass Specific Metabolic Rate") + ylab("Deviance Residuals")
        
    QQNorm.1.m1 = ggplot(data.frame("Res" = residuals(AAS.MS.m1), "Fit" = predict(AAS.MS.m1)),
        aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone.m1, Res.by.Fit.m1, Res.by.Resp.m1, QQNorm.1.m1, nrow = 2, top = 
        "Mass Specific Resting Metabolic Rate Residuals")

## Awfully course, but fair residuals. Double-checking normality.

    qqPlot(residuals(AAS.MS.m1, type = "deviance"), ylab = "Deviance Residuals")

## 172  83 
## 129  63
    # Reasonable

r.squaredGLMM(AAS.MS.m1)
##             R2m        R2c
## [1,] 0.06030224 0.06030224
# As expected. Remarkably low.

# No clear rationale to test variance structure. Calculating p-values.
 
AAS.dwdata1 = data.frame(term = c(rownames(summary(AAS.MS.m1)$coefficients)),
                         estimate = c(summary(AAS.MS.m1)$coefficients[1:2]),
                         std.error = c(summary(AAS.MS.m1)$coefficients[3:4]),
                         statistic = c(summary(AAS.MS.m1)$coefficients[5:6]))

p.values = c()
for (i in 1:nrow(AAS.dwdata1)){
  p.values[i] = pt(abs(AAS.dwdata1$statistic[i]), df = 8, lower.tail = FALSE)
} 
AAS.dwdata1$p.values = round(p.values, 3)

AAS.dwdata1$term = c("(Intercept)", 
                     "Offspring Acclimation Temperature")

print(AAS.dwdata1)
##                                term   estimate  std.error  statistic p.values
## 1                       (Intercept)  5.5645121 0.03302392 168.499417    0.000
## 2 Offspring Acclimation Temperature -0.1494679 0.04538704  -3.293184    0.005

###Testing Model Assumptions for Absolute Aerobic Scope Models: Model 2.

AAS.MS.m2 = lmer(log(AAS.MS) ~ Mother.Acc + Offspring.Acc + (1|Mother.ID) + (1|Father.ID), 
        REML = FALSE, data = AAS.MS.Data)

summary(AAS.MS.m2)$varcor 
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.0000  
##  Father.ID (Intercept) 0.0000  
##  Residual              0.2949
# Again, Father ID and Mother ID alone explained zero residual variance.

    Res.Alone.m2 = ggplot(data = data.frame("Res" = residuals(AAS.MS.m2, type = "deviance")), 
        aes(x = 1:nrow(AAS.MS.Data), y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
        ylab("Deviance Residuals")

    Res.by.Fit.m2 = ggplot(data = data.frame("Res" = residuals(AAS.MS.m2), "Fit" = predict(AAS.MS.m2)), 
        aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
        ylab("Deviance Residuals")

    Res.by.Resp.m2 = ggplot(data = data.frame("Res" = residuals(AAS.MS.m2), "Resp" = log(AAS.MS.Data$AAS.MS)), 
        aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
        xlab("log Mass Specific Metabolic Rate") + ylab("Deviance Residuals")
        
    QQNorm.1.m2 = ggplot(data.frame("Res" = residuals(AAS.MS.m2), "Fit" = predict(AAS.MS.m2)),
        aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone.m2, Res.by.Fit.m2, Res.by.Resp.m2, QQNorm.1.m2, nrow = 2, top = 
        "Mass Specific Resting Metabolic Rate Residuals")

## Residuals similar to model 1, with slight right tail.

    qqPlot(residuals(AAS.MS.m2, type = "deviance"), ylab = "Deviance Residuals")

## 172  83 
## 129  63
# As above.

# Parameter-specific heteroskedasticity

    {
    DamT.Res = ggplot(data = AAS.MS.Data, aes(x = as.factor(Mother.Acc), 
        y = residuals(AAS.MS.m2, "deviance"))) + my.theme + theme_bw() + 
        geom_boxplot(fill = "mediumorchid", aes(middle = mean(residuals(AAS.MS.m2, "deviance")))) + 
        xlab("Dam Acclimation Temperature") + ylab("Deviance Residuals")

    OffT.Res = ggplot(data = AAS.MS.Data, aes(x = as.factor(Offspring.Acc), 
        y = residuals(AAS.MS.m2, "deviance"))) + my.theme + theme_bw() + 
        geom_boxplot(fill = "mediumseagreen", aes(middle = mean(residuals(AAS.MS.m2, "deviance")))) + 
        xlab("Offspring Acclimation Temperature") + ylab("Deviance Residuals")

    Off.by.Dam.T.Res = ggplot(data = AAS.MS.Data, aes(x = as.factor(Offspring.Acc), 
        y = residuals(AAS.MS.m2, type = "pearson"), fill = Mother.Acc)) + 
        my.theme + theme_bw() + geom_boxplot() + scale_fill_manual(values = c("mediumorchid", 
        "mediumseagreen")) + xlab("Offspring Acclimation Temperature") + ylab("Deviance Residuals") + 
      theme(legend.position = "top")
    }

    grid.arrange(DamT.Res, OffT.Res, Off.by.Dam.T.Res, nrow = 1, 
    top = "Fixed Effects by Residuals: Boxplot Bars Represent Means")  

    # Probable offspring by maternal heteroskedasticity.

    AAS.MS.m2.2 = lme(log(AAS.MS) ~ Mother.Acc + Offspring.Acc, random = list(~1|Mother.ID, ~1|Father.ID),
        weights = varIdent(form = ~Offspring.Acc|Mother.Acc), data = AAS.MS.Data)

    lapply(list(AAS.MS.m2, AAS.MS.m2.2), r.squaredGLMM)
## [[1]]
##             R2m        R2c
## [1,] 0.06331924 0.06331924
## 
## [[2]]
##             R2m        R2c
## [1,] 0.06904318 0.06904318
    # Negligable variance gained. Maintaining first model.

# Calculating p-values

AAS.dwdata2 = data.frame(term = c(rownames(summary(AAS.MS.m2)$coefficients)),
                         estimate = c(summary(AAS.MS.m2)$coefficients[1:3]),
                         std.error = c(summary(AAS.MS.m2)$coefficients[4:6]),
                         statistic = c(summary(AAS.MS.m2)$coefficients[7:9]))

p.values = c()
for (i in 1:nrow(AAS.dwdata2)){
  p.values[i] = pt(abs(AAS.dwdata2$statistic[i]), df = 8, lower.tail = FALSE)
} 
AAS.dwdata2$p.values = round(p.values, 3)

AAS.dwdata2$term = c("(Intercept)", "Dam Acclimation Temperature", 
                     "Offspring Acclimation Temperature")

print(AAS.dwdata2) 
##                                term    estimate  std.error   statistic p.values
## 1                       (Intercept)  5.54564357 0.04172416 132.9120402    0.000
## 2       Dam Acclimation Temperature  0.03354401 0.04545720   0.7379251    0.241
## 3 Offspring Acclimation Temperature -0.14848951 0.04533392  -3.2754615    0.006

###Testing Model Assumptions for Absolute Aerobic Scope Models: Model 3.

AAS.MS.m3 = lmer(log(AAS.MS) ~ Father.Acc + Offspring.Acc + 
                   (1|Mother.ID) + (1|Father.ID), REML = FALSE, data = AAS.MS.Data)

summary(AAS.MS.m3)$varcor 
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.00000 
##  Father.ID (Intercept) 0.00000 
##  Residual              0.29495
# Father ID alone explained zero variance.

# Residual plots
    Res.Alone.m2 = ggplot(data = data.frame("Res" = residuals(AAS.MS.m3, type = "deviance")), 
        aes(x = 1:nrow(AAS.MS.Data), y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
        ylab("Deviance Residuals")

    Res.by.Fit.m2 = ggplot(data = data.frame("Res" = residuals(AAS.MS.m3), "Fit" = predict(AAS.MS.m3)), 
        aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
        ylab("Deviance Residuals")

    Res.by.Resp.m2 = ggplot(data = data.frame("Res" = residuals(AAS.MS.m3), "Resp" = log(AAS.MS.Data$AAS.MS)), 
        aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
        xlab("log Mass Specific Metabolic Rate") + ylab("Deviance Residuals")
        
    QQNorm.1.m2 = ggplot(data.frame("Res" = residuals(AAS.MS.m3), "Fit" = predict(AAS.MS.m3)),
        aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone.m3, Res.by.Fit.m3, Res.by.Resp.m3, QQNorm.1.m3, nrow = 2, top = 
        "Mass Specific Resting Metabolic Rate Residuals")

## Residuals straying.

    qqPlot(residuals(AAS.MS.m3, type = "deviance"), ylab = "Deviance Residuals")

## 172  83 
## 129  63
    # Okay, still within range.

r.squaredGLMM(AAS.MS.m3)
##             R2m        R2c
## [1,] 0.06304827 0.06304827
# As low as previous two models.

    {
    PatT.Res = ggplot(data = AAS.MS.Data, aes(x = as.factor(Father.Acc), 
        y = residuals(AAS.MS.m2, "deviance"))) + my.theme + theme_bw() + 
        geom_boxplot(fill = "mediumorchid", aes(middle = mean(residuals(AAS.MS.m2, "deviance")))) + 
        xlab("Dam Acclimation Temperature") + ylab("Deviance Residuals")

    OffT.Res = ggplot(data = AAS.MS.Data, aes(x = as.factor(Offspring.Acc), 
        y = residuals(AAS.MS.m2, "deviance"))) + my.theme + theme_bw() + 
        geom_boxplot(fill = "mediumseagreen", aes(middle = mean(residuals(AAS.MS.m2, "deviance")))) + 
        xlab("Offspring Acclimation Temperature") + ylab("Deviance Residuals")

    Off.by.Pat.T.Res = ggplot(data = AAS.MS.Data, aes(x = as.factor(Offspring.Acc), 
        y = residuals(AAS.MS.m2.2, type = "pearson"), fill = Father.Acc)) + 
        my.theme + theme_bw() + geom_boxplot() + scale_fill_manual(values = c("mediumorchid", 
        "mediumseagreen")) + xlab("Offspring Acclimation Temperature") + ylab("Deviance Residuals") + 
      theme(legend.position = "top")
    }
    grid.arrange(PatT.Res, OffT.Res, Off.by.Pat.T.Res, nrow = 1, 
    top = "Fixed Effects by Residuals: Boxplot Bars Represent Means") 

## No cause for concern.

# Model result plots

AAS.dwdata3 = data.frame(term = c(rownames(summary(AAS.MS.m3)$coefficients)),
                         estimate = c(summary(AAS.MS.m3)$coefficients[1:3]),
                         std.error = c(summary(AAS.MS.m3)$coefficients[4:6]),
                         statistic = c(summary(AAS.MS.m3)$coefficients[7:9]))

p.values = c()
for (i in 1:nrow(AAS.dwdata3)){
  p.values[i] = pt(abs(AAS.dwdata3$statistic[i]), df = 8, lower.tail = FALSE)
} 
AAS.dwdata3$p.values = round(p.values, 3)

AAS.dwdata3$term = c("(Intercept)", "Sire Acclimation Temperature",
                     "Offspring Acclimation Temperature")

## Little reason to overlap model outputs, given small number of fixed effects in each.

AAS.dwdata1 = AAS.dwdata1 %>% mutate(model = "Model 1")
AAS.dwdata2 = AAS.dwdata2 %>% mutate(model = "Model 2")
AAS.dwdata3 = AAS.dwdata3 %>% mutate(model = "Model 3")
AAS.MS.Models = rbind(AAS.dwdata1, AAS.dwdata2, AAS.dwdata3)

AAS.Fit.m1 = data.frame("Response" = c(log(AAS.MS.Data$AAS.MS)), "Fit" = c(predict(AAS.MS.m1)))
AAS.Fit.m2 = data.frame("Response" = c(log(AAS.MS.Data$AAS.MS)), "Fit" = c(predict(AAS.MS.m2)))
AAS.Fit.m3 = data.frame("Response" = c(log(AAS.MS.Data$AAS.MS)), "Fit" = c(predict(AAS.MS.m3)))
AAS.Fit.m1 = AAS.Fit.m1 %>% mutate(model = "Model 1")
AAS.Fit.m2 = AAS.Fit.m2 %>% mutate(model = "Model 2")
AAS.Fit.m3 = AAS.Fit.m3 %>% mutate(model = "Model 3")

AAS.Fit.All = rbind(AAS.Fit.m1, AAS.Fit.m2, AAS.Fit.m3)
AAS.Fits = ggplot(data = AAS.Fit.All, aes(x = Fit, y = Response, fill = model, colour = model,
           linetype = model)) + theme_bw() + my.theme + geom_point(pch=21, size = 3, alpha = 0.5) + 
           geom_smooth(method = "lm") + xlab("Predicted log Absolute Aerobic Scope") + 
           ylab("Measured log Absolute \n Aerobic Scope") + scale_fill_manual(values = c("mediumseagreen", 
           "mediumorchid", "cornflowerblue")) + scale_colour_manual(values = c("black", 
         "black", "black")) + theme(legend.title=element_blank())

print(AAS.Fits)

ggsave("/home/joshk/git_repositories/BrookTrout_TG/Figures/Absolute Aerobic Scope Model Fit Comparisons.jpeg", dpi = 800, limitsize = TRUE)

##Non-Mass Specific Model Assessment ###Testing Model Assumptions for Criticial Thermal Maximum: Model 1.

## First assessing response by predictors;

    {
    PatT.CTM = ggplot(data = CTM.NMS.Data, aes(x = as.factor(Father.Acc), 
        y = CTM)) + theme_bw() + 
        geom_boxplot(fill = "mediumorchid", aes(middle = mean(CTM))) + 
        xlab("Sire Acclimation Temperature") + ylab("Critical Thermal Maximum \n (Degrees C)")

    DamT.CTM  = ggplot(data = CTM.NMS.Data, aes(x = as.factor(Mother.Acc), 
        y = CTM)) + theme_bw() + 
        geom_boxplot(fill = "mediumseagreen", aes(middle = mean(CTM))) + 
        xlab("Dam Acclimation Temperature") + ylab("Critical Thermal Maximum \n (Degrees C)")

    OffT.CTM  = ggplot(data = CTM.NMS.Data, aes(x = as.factor(Offspring.Acc), 
        y = CTM)) + theme_bw() + 
        geom_boxplot(fill = "royalblue", aes(middle = mean(CTM))) + 
        xlab("Offspring Acclimation Temperature") + ylab("Critical Thermal Maximum \n (Degrees C)")

    Pat.by.CTM1  = ggplot(data = CTM.NMS.Data, aes(x = as.factor(Father.Acc), 
        y = CTM, fill = Offspring.Acc)) + theme_bw() + 
        theme(legend.position = "top") + geom_boxplot() + 
        scale_fill_manual(values = c("mediumorchid", "royalblue")) + 
        xlab("Sire Acclimation Temperature") + ylab("Critical Thermal Maximum \n (Degrees C)")

    Pat.by.CTM2 = ggplot(data = CTM.NMS.Data, aes(x = as.factor(Father.Acc), 
        y = CTM, fill = Mother.Acc)) + theme_bw() + 
        theme(legend.position = "top") + geom_boxplot() + scale_fill_manual(values = c("mediumseagreen", 
        "royalblue")) + xlab("Sire Acclimation Temperature") + ylab("Critical Thermal Maximum \n (Degrees C)")

    Off.by.CTM = ggplot(data = CTM.NMS.Data, aes(x = as.factor(Offspring.Acc), 
        y = CTM, fill = Mother.Acc)) + theme_bw() + 
        theme(legend.position = "top") + geom_boxplot() + scale_fill_manual(values = c("mediumorchid", 
        "mediumseagreen")) + xlab("Offspring Acclimation Temperature") + 
      ylab("Critical Thermal Maximum \n (Degrees C)")

    }

    grid.arrange(PatT.CTM, DamT.CTM, OffT.CTM, Pat.by.CTM1, Pat.by.CTM2,
    Off.by.CTM, nrow = 2, 
    top = "Fixed Effects by Response: Boxplot Bars Represent Means") 

## MODEL

CTM.m1 = lmer(CTM ~ Father.Acc + Mother.Acc + Offspring.Acc + Father.Acc:Mother.Acc + 
        Father.Acc:Offspring.Acc + Mother.Acc:Offspring.Acc + 
        Father.Acc:Mother.Acc:Offspring.Acc + (1|Mother.ID) + (1|Father.ID), 
        REML = FALSE, weights = Mass, data = CTM.NMS.Data)

summary(CTM.m1)$varcor 
##  Groups    Name        Std.Dev.  
##  Mother.ID (Intercept) 2.0348e-02
##  Father.ID (Intercept) 9.7330e-09
##  Residual              4.0160e-01
## Paternal ID dropping from model.

    Res.Alone.m1 = ggplot(data = data.frame("Res" = residuals(CTM.m1, type = "deviance")), 
        aes(x = 1:nrow(CTM.NMS.Data), y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
        ylab("Deviance Residuals")

    Res.by.Fit.m1 = ggplot(data = data.frame("Res" = residuals(CTM.m1), "Fit" = predict(CTM.m1)), 
        aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
        ylab("Deviance Residuals")

    Res.by.Resp.m1 = ggplot(data = data.frame("Res" = residuals(CTM.m1), "Resp" = CTM.NMS.Data$CTM), 
        aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
        geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
        xlab("log Mass Specific Metabolic Rate") + ylab("Deviance Residuals")
        
    QQNorm.1.m1 = ggplot(data.frame("Res" = residuals(CTM.m1), "Fit" = predict(CTM.m1)),
        aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone.m1, Res.by.Fit.m1, Res.by.Resp.m1, QQNorm.1.m1, nrow = 2, top = 
        "Mass Specific Resting Metabolic Rate Residuals")

    qqPlot(residuals(CTM.m1, type = "deviance"), ylab = "Deviance Residuals")

## 102  14 
##  97  14
    # Some low end deviation.
            
CTM.m1.res = data.frame("Model" = "Model1", "Res" = c(residuals(CTM.m1, type = "deviance")))

ggplot(CTM.m1.res, aes(Res, fill = Model, colour = Model)) + geom_histogram(alpha = 0.5) + theme_bw() + 
    my.theme + xlab("Deviance Residuals") + ylab("Count") + scale_fill_manual(values = 
    c("mediumseagreen")) + scale_colour_manual(values = c("black")) + theme(legend.position = "none") + 
    geom_density(alpha = 0.5, aes(y=0.1 * ..count..))

# Reasonable! 
    
## Assessing parameter-specific patterns.

    {
    PatT.Res.m1 = ggplot(data = CTM.NMS.Data, aes(x = as.factor(Father.Acc), 
        y = residuals(CTM.m1, "deviance"))) + theme_bw() + 
        geom_boxplot(fill = "mediumorchid", aes(middle = mean(residuals(CTM.m1, "deviance")))) + 
        xlab("Sire Acclimation Temperature") + ylab("Deviance Residuals")

    DamT.Res.m1 = ggplot(data = CTM.NMS.Data, aes(x = as.factor(Mother.Acc), 
        y = residuals(CTM.m1, "deviance"))) + theme_bw() + 
        geom_boxplot(fill = "mediumseagreen", aes(middle = mean(residuals(CTM.m1, "deviance")))) + 
        xlab("Dam Acclimation Temperature") + ylab("Deviance Residuals")

    OffT.Res.m1 = ggplot(data = CTM.NMS.Data, aes(x = as.factor(Offspring.Acc), 
        y = residuals(CTM.m1, "deviance"))) + theme_bw() + 
        geom_boxplot(fill = "royalblue", aes(middle = mean(residuals(CTM.m1, "deviance")))) + 
        xlab("Offspring Acclimation Temperature") + ylab("Deviance Residuals")

    Pat.by.Off.T.Res.m1 = ggplot(data = CTM.NMS.Data, aes(x = as.factor(Father.Acc), 
        y = residuals(CTM.m1, type = "deviance"), fill = Offspring.Acc)) + theme_bw() + 
        theme(legend.position = "top") + geom_boxplot() + 
        scale_fill_manual(values = c("mediumorchid", "royalblue")) + 
        xlab("Sire Acclimation Temperature") + ylab("Deviance Residuals")

    Pat.by.Dam.T.Res.m1 = ggplot(data = CTM.NMS.Data, aes(x = as.factor(Father.Acc), 
        y = residuals(CTM.m1, type = "deviance"), fill = Mother.Acc)) + theme_bw() + 
        theme(legend.position = "top") + geom_boxplot() + scale_fill_manual(values = c("mediumseagreen", 
        "royalblue")) + xlab("Sire Acclimation Temperature") + ylab("Deviance Residuals")

    Off.by.Dam.T.Res.m1 = ggplot(data = CTM.NMS.Data, aes(x = as.factor(Offspring.Acc), 
        y = residuals(CTM.m1, type = "deviance"), fill = Mother.Acc)) + theme_bw() + 
        theme(legend.position = "top") + geom_boxplot() + scale_fill_manual(values = c("mediumorchid", 
        "mediumseagreen")) + xlab("Offspring Acclimation Temperature") + ylab("Deviance Residuals")

    }
    grid.arrange(PatT.Res.m1, DamT.Res.m1, OffT.Res.m1, Pat.by.Off.T.Res.m1, Pat.by.Dam.T.Res.m1,
    Off.by.Dam.T.Res.m1, nrow = 2, 
    top = "Fixed Effects by Residuals: Boxplot Bars Represent Means") 

# Possible sire by offspring heteroskedasticity. Retrying model with constant variance per group.

    CTM.m1.1 = lme(CTM ~ Father.Acc + Mother.Acc + Offspring.Acc + Father.Acc:Mother.Acc + 
        Father.Acc:Offspring.Acc + Mother.Acc:Offspring.Acc + 
        Father.Acc:Mother.Acc:Offspring.Acc, random = list(~1|Mother.ID, ~1|Father.ID), 
        weights = ~Mass, data = CTM.NMS.Data)

    CTM.m1.2 = lme(CTM ~ Father.Acc + Mother.Acc + Offspring.Acc + Father.Acc:Mother.Acc + 
        Father.Acc:Offspring.Acc + Mother.Acc:Offspring.Acc + 
        Father.Acc:Mother.Acc:Offspring.Acc, random = list(~1|Mother.ID, ~1|Father.ID), 
        weights = varIdent(form = ~Mass|Father.Acc*Offspring.Acc), data = CTM.NMS.Data)

summary(CTM.m1.1)$AIC/summary(CTM.m1.2)$AIC
## [1] 6.157988
    # Wow! Significant improvement in AIC. Residuals?

    Res.Alone.m1 = ggplot(data = data.frame("Res" = residuals(CTM.m1.2, type = "normalized")), 
        aes(x = 1:nrow(CTM.NMS.Data), y = Res)) + theme_bw() + 
        geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
        ylab("Deviance Residuals")

    Res.by.Fit.m1 = ggplot(data = data.frame("Res" = residuals(CTM.m1.2, type = "normalized"), 
        "Fit" = predict(CTM.m1.2)), aes(x = Fit, y = Res)) + theme_bw() + 
        geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
        ylab("Deviance Residuals")

    Res.by.Resp.m1 = ggplot(data = data.frame("Res" = residuals(CTM.m1.2, type = "normalized"), 
        "Resp" = CTM.NMS.Data$CTM), aes(x = Resp, y = Res)) + theme_bw() + 
        geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
        xlab("log Mass Specific Metabolic Rate") + ylab("Deviance Residuals")
        
    QQNorm.1.m1 = ggplot(data.frame("Res" = residuals(CTM.m1.2, type = "normalized"), 
        "Fit" = predict(CTM.m1.2)), aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + 
         theme_bw()

    grid.arrange(Res.Alone.m1, Res.by.Fit.m1, Res.by.Resp.m1, QQNorm.1.m1, nrow = 2, top = 
        "Mass Specific Resting Metabolic Rate Residuals")

    # Far cleaner. Double-checking confidence intervals.

    qqPlot(residuals(CTM.m1.2, type = "normalized"), ylab = "Normalized Residuals")

## 1/3 1/3 
## 154 161
    # Well enough behaved! Re-testing AIC with new variance structure.

    CTM.model = lme(CTM ~ Mother.Acc*Father.Acc*Offspring.Acc + Mass, random = list(~1|Mother.ID, ~1|Father.ID), 
    weights = varIdent(~Mass|Father.Acc*Offspring.Acc), data = CTM.NMS.Data)
    
    subset(dredge(CTM.model, rank = "AICc"), delta <= 5.0)
## Global model call: lme.formula(fixed = CTM ~ Mother.Acc * Father.Acc * Offspring.Acc + 
##     Mass, data = CTM.NMS.Data, random = list(~1 | Mother.ID, 
##     ~1 | Father.ID), weights = varIdent(~Mass | Father.Acc * 
##     Offspring.Acc))
## ---
## Model selection table 
##    (Int) Mth.Acc Off.Acc Mth.Acc:Off.Acc df  logLik AICc delta weight
## 9  28.57               +                  5 -19.117 48.5  0.00   0.79
## 77 28.48       +       +               +  7 -18.799 52.1  3.62   0.13
## 13 28.54       +       +                  6 -20.346 53.1  4.57   0.08
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Mother.ID', '1 | Father.ID %in% Mother.ID'

##Re-calculating AICs for Model Iterations: Criticial Thermal Maximum.

AIC.Results.CTM1 = subset(dredge(CTM.model, rank = "AICc", 
    extra = list(R2 = function(x) r.squaredGLMM(x))), delta < 5)
Select.Models.CTM = c()
Delta.NMS = c()
Response = c()
Marg.R2 = c() 
Cond.R2 = c()
Select.Models.CTM1 = gsub("[[:space:]].*", "", gsub(".*~ ", "", paste(attr(AIC.Results.CTM1, "model.calls")[[1]])[2]))
Delta.CTM1 = AIC.Results.CTM1$delta[1]
Response = gsub(" ~.*", "", paste(attr(AIC.Results.CTM1, "model.calls")[[1]])[2])
Marg.R2 = AIC.Results.CTM1$R21[1]
Cond.R2 = AIC.Results.CTM1$R22[1]
Table = data.frame("Reponse.Type" = "Mass Independant", "Response" = Response, 
"Delta.AIC" = Delta.CTM1, "Marginal.R2" = Marg.R2, "Conditional.R2" = Cond.R2,
"Equation" = Select.Models.CTM1)
assign(paste("AIC.Selection", "NMS", 4, sep = "_"), Table)

# Stacking all top models

AIC.Selection.NMS = rbind(AIC.Selection_NMS_1, AIC.Selection_NMS_2, AIC.Selection_NMS_3, AIC.Selection_NMS_4)
AIC.Selection = rbind(AIC.Selection.MS, AIC.Selection.NMS)
print(AIC.Selection, right = F)
##    Reponse.Type     Response    Delta.AIC Marginal.R2 Conditional.R2
## 1  Mass-Specific    log(RMR.MS) 0.0000000 0.09689341  0.23258538    
## 2  Mass-Specific    log(RMR.MS) 0.4320103 0.05317331  0.24043771    
## 3  Mass-Specific    log(RMR.MS) 1.2633111 0.11949983  0.22942058    
## 4  Mass-Specific    log(RMR.MS) 1.6197476 0.07652765  0.23332101    
## 5  Mass-Specific    log(RMR.MS) 2.1592505 0.09688383  0.23257954    
## 6  Mass-Specific    log(RMR.MS) 3.4158810 0.11948859  0.22852419    
## 7  Mass-Specific    log(RMR.MS) 3.4339163 0.11966268  0.23078182    
## 8  Mass-Specific    log(RMR.MS) 3.4465742 0.11949854  0.22941555    
## 9  Mass-Specific    log(RMR.MS) 3.7625995 0.07663967  0.23491827    
## 10 Mass-Specific    log(MMR.MS) 0.0000000 0.12081940  0.14201563    
## 11 Mass-Specific    log(MMR.MS) 1.1337870 0.12482063  0.14707607    
## 12 Mass-Specific    log(MMR.MS) 1.6198084 0.07304761  0.14942805    
## 13 Mass-Specific    log(MMR.MS) 1.6632302 0.12283706  0.14581838    
## 14 Mass-Specific    log(MMR.MS) 2.7260488 0.07701051  0.15563521    
## 15 Mass-Specific    log(MMR.MS) 2.7981663 0.12640953  0.14714116    
## 16 Mass-Specific    log(MMR.MS) 2.8091275 0.12687437  0.15094298    
## 17 Mass-Specific    log(MMR.MS) 3.1451250 0.12529971  0.14712299    
## 18 Mass-Specific    log(MMR.MS) 4.4679396 0.07829618  0.15562601    
## 19 Mass-Specific    log(MMR.MS) 4.5445400 0.12823397  0.15073801    
## 20 Mass-Specific    log(MMR.MS) 4.7897724 0.12703622  0.14722626    
## 21 Mass-Specific    log(MMR.MS) 4.8318359 0.12737623  0.15100909    
## 22 Mass-Specific    log(AAS.MS) 0.0000000 0.06030224  0.06030224    
## 23 Mass-Specific    log(AAS.MS) 1.6058206 0.06331924  0.06331924    
## 24 Mass-Specific    log(AAS.MS) 1.6547206 0.06304827  0.06304827    
## 25 Mass-Specific    log(AAS.MS) 2.7150419 0.06921105  0.06921105    
## 26 Mass-Specific    log(AAS.MS) 3.2960644 0.06600676  0.06600676    
## 27 Mass-Specific    log(AAS.MS) 3.8171171 0.06312377  0.06312377    
## 28 Mass-Specific    log(AAS.MS) 4.4085037 0.07201229  0.07201229    
## 29 Mass Independant log(RMR)    0.0000000 0.13398490  0.18729777    
## 30 Mass Independant log(RMR)    1.2602254 0.14177631  0.20537453    
## 31 Mass Independant log(RMR)    1.6073708 0.15042976  0.21764843    
## 32 Mass Independant log(RMR)    1.6474074 0.12869909  0.17308051    
## 33 Mass Independant log(RMR)    2.7321064 0.13274864  0.17071098    
## 34 Mass Independant log(RMR)    2.8878776 0.13578794  0.19007496    
## 35 Mass Independant log(RMR)    3.2743628 0.14463523  0.20289714    
## 36 Mass Independant log(RMR)    4.0588194 0.14000948  0.18880081    
## 37 Mass Independant log(RMR)    4.0950899 0.09988873  0.21612058    
## 38 Mass Independant log(RMR)    4.5640317 0.14806374  0.20017371    
## 39 Mass Independant log(RMR)    4.8488012 0.13927880  0.19821267    
## 40 Mass Independant log(MMR)    0.0000000 0.36920443  0.37186427    
## 41 Mass Independant log(MMR)    0.5548350 0.37942598  0.38649376    
## 42 Mass Independant log(MMR)    0.6787176 0.36963415  0.36963415    
## 43 Mass Independant log(MMR)    0.7216353 0.37510675  0.37928176    
## 44 Mass Independant log(MMR)    0.8325546 0.37557342  0.37704611    
## 45 Mass Independant log(MMR)    1.2075921 0.38568092  0.39455042    
## 46 Mass Independant log(MMR)    1.3053415 0.37389091  0.37398454    
## 47 Mass Independant log(MMR)    1.4075558 0.38019832  0.38260583    
## 48 Mass Independant log(MMR)    2.6399397 0.37676558  0.37901010    
## 49 Mass Independant log(MMR)    3.2056969 0.38148663  0.38472604    
## 50 Mass Independant log(MMR)    3.2512482 0.37444184  0.37444184    
## 51 Mass Independant log(MMR)    3.3253542 0.38543858  0.39392723    
## 52 Mass Independant log(MMR)    3.3308469 0.38082566  0.38300719    
## 53 Mass Independant log(MMR)    3.5169143 0.38027035  0.38224935    
## 54 Mass Independant AAS         0.0000000 0.37930675  0.37930675    
## 55 Mass Independant AAS         1.5094369 0.38174466  0.38174466    
## 56 Mass Independant AAS         2.1236991 0.37949846  0.37949846    
## 57 Mass Independant AAS         2.6531066 0.38559940  0.38559940    
## 58 Mass Independant AAS         3.6519021 0.38196577  0.38196577    
## 59 Mass Independant AAS         4.3197068 0.37952426  0.37952426    
## 60 Mass Independant AAS         4.8437724 0.38574420  0.38574420    
## 61 Mass Independant CTM         0.0000000 0.46835156  0.50879302    
##    Equation                                                                                                                                  
## 1  Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 2  Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                         
## 3  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                               
## 4  Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 5  Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                 
## 6  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                       
## 7  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                    
## 8  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                    
## 9  Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                 
## 10 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 11 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                               
## 12 Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                         
## 13 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                 
## 14 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 15 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                    
## 16 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                    
## 17 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                       
## 18 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                 
## 19 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc + Mother.Acc:Offspring.Acc         
## 20 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc + Mother.Acc:Offspring.Acc            
## 21 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc + Father.Acc:Offspring.Acc            
## 22 Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                         
## 23 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 24 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 25 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                 
## 26 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                               
## 27 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                 
## 28 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                    
## 29 log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                
## 30 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 31 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                     
## 32 Father.Acc + log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                   
## 33 Father.Acc + log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                           
## 34 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                   
## 35 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                        
## 36 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                           
## 37 log(Mass) + (1 | Mother.ID) + (1 | Father.ID)                                                                                             
## 38 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc + Mother.Acc:Offspring.Acc
## 39 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                        
## 40 log(Mass) + (1 | Mother.ID) + (1 | Father.ID)                                                                                             
## 41 log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                             
## 42 Father.Acc + log(Mass) + (1 | Mother.ID) + (1 | Father.ID)                                                                                
## 43 log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                
## 44 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 45 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 46 Father.Acc + log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                   
## 47 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                   
## 48 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                     
## 49 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                        
## 50 Father.Acc + log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                           
## 51 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                     
## 52 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                           
## 53 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                        
## 54 log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                             
## 55 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 56 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 57 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                     
## 58 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                   
## 59 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                     
## 60 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                        
## 61 Offspring.Acc
write.csv(AIC.Selection, "Model_Outputs_AICc.csv", row.names = FALSE)

Reassessing Criticial Thermal Maximum: Model 1.

    CTM.model.m1 = lme(CTM ~ Offspring.Acc, random = list(~1|Mother.ID, ~1|Father.ID), 
    weights = varIdent(~Mass|Father.Acc*Offspring.Acc), data = CTM.NMS.Data)

    # Residual Plots

    Res.Alone.m1 = ggplot(data = data.frame("Res" = residuals(CTM.model.m1, type = "normalized")), 
        aes(x = 1:nrow(CTM.NMS.Data), y = Res)) + theme_bw() + 
        geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
        ylab("Deviance Residuals")

    Res.by.Fit.m1 = ggplot(data = data.frame("Res" = residuals(CTM.model.m1, type = "normalized"), 
        "Fit" = predict(CTM.model.m1)), aes(x = Fit, y = Res)) + theme_bw() + 
        geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
        ylab("Deviance Residuals")

    Res.by.Resp.m1 = ggplot(data = data.frame("Res" = residuals(CTM.model.m1, type = "normalized"), 
        "Resp" = CTM.NMS.Data$CTM), aes(x = Resp, y = Res)) + theme_bw() + 
        geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
        xlab("log Mass Specific Metabolic Rate") + ylab("Deviance Residuals")
        
    QQNorm.1.m1 = ggplot(data.frame("Res" = residuals(CTM.model.m1, type = "normalized"), 
        "Fit" = predict(CTM.model.m1)), aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + 
         theme_bw()

    grid.arrange(Res.Alone.m1, Res.by.Fit.m1, Res.by.Resp.m1, QQNorm.1.m1, nrow = 2, top = 
        "Mass Specific Resting Metabolic Rate Residuals")

    qqPlot(residuals(CTM.model.m1, type = "normalized"), ylab = "Deviance Residuals")

## 1/1 3/3 
##  14  89
    r.squaredGLMM(CTM.model.m1)
##            R2m      R2c
## [1,] 0.4683516 0.508793
    # Surprisingly reasonable.

# Calculating P-values

CTM.dwdata1 = data.frame(term = c(names(summary(CTM.model.m1)$coefficients$fixed)),
                         estimate = c(summary(CTM.model.m1)$coefficients$fixed[1:2]),
                         std.error = c(summary(CTM.model.m1)$tTable[,2]),
                         statistic = c(summary(CTM.model.m1)$tTable[,4]))

p.values = c()
for (i in 1:nrow(CTM.dwdata1)){
  p.values[i] = pt(abs(CTM.dwdata1$statistic[i]), df = 8, lower.tail = FALSE)
} 

CTM.dwdata1$p.values = round(p.values, 3)

CTM.dwdata1$term = c("(Intercept)", "Offspring Acclimation Temperature")

# Plotting fitted model

CTM.Fits = ggplot(data = data.frame(Pred = c(predict(CTM.model.m1)), Resp = c(CTM.NMS.Data$CTM), 
                                    Model = "Mod1"), 
        aes(x = Pred, y = Resp, fill = Model, colour = Model)) + theme_bw() + my.theme + 
        geom_point(pch=21, size = 3, alpha = 0.5) + geom_smooth(method = "lm") + 
        xlab("Predicted Critical Thermal Maximum") + ylab("Measured Critical \n Thermal Maximum") + 
        scale_fill_manual(values = c("cornflowerblue")) + scale_colour_manual(values = c("black")) + 
        theme(legend.position = "none") + xlim(c(28.4,29.2))

print(CTM.Fits)

ggsave("/home/joshk/git_repositories/BrookTrout_TG/Figures/Critical Thermal Maximum Fit Comparisons.jpeg", dpi = 800, limitsize = TRUE)

CTM.Models = CTM.dwdata1 %>% mutate(model = "Model 1")

##Almagamating Above Results.

require('plyr')

R2.Dataframe = ldply(lapply(list(RMR.MS.m1, RMR.MS.m2, RMR.MS.m3, 
                                 MMR.MS.m1, MMR.MS.m2, MMR.MS.m3, 
                                 AAS.MS.m1, AAS.MS.m2, AAS.MS.m3,
                                 CTM.m1), r.squaredGLMM), data.frame)

RMRDF1 = subset(RMR.MS.Models, model == "Model 1") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 1")%>% 
  mutate("Marginal R2" = R2.Dataframe[1,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[1,2])%>%
  mutate("AIC" = summary(RMR.MS.m1)$AICtab[1])
RMRDF2 = subset(RMR.MS.Models, model == "Model 2") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 2")%>% 
  mutate("Marginal R2" = R2.Dataframe[2,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[2,2]) %>%
  mutate("AIC" = summary(RMR.MS.m2)$AICtab[1]) 
RMRDF3 = subset(RMR.MS.Models, model == "Model 3") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 3")%>% 
  mutate("Marginal R2" = R2.Dataframe[3,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[3,2])%>% 
  mutate("AIC" = summary(RMR.MS.m3)$AICtab[1])

RMRDF = rbind(RMRDF1, RMRDF2, RMRDF3)

RMRDF$Response = "log Resting Metabolic Rate" 
RMRDF$Response = "log Resting Metabolic Rate" 
RMRDF$Response = "log Resting Metabolic Rate" 

MMRDF1 = subset(MMR.MS.Models, model == "Model 1") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 1")%>% 
  mutate("Marginal R2" = R2.Dataframe[4,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[4,2]) %>% 
  mutate("AIC" = summary(MMR.MS.m1)$AICtab[1])
MMRDF2 = subset(MMR.MS.Models, model == "Model 2") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 2")%>% 
  mutate("Marginal R2" = R2.Dataframe[5,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[5,2]) %>% 
  mutate("AIC" = summary(MMR.MS.m2)$AICtab[1])
MMRDF3 = subset(MMR.MS.Models, model == "Model 3") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 3")%>% 
  mutate("Marginal R2" = R2.Dataframe[6,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[6,2]) %>% 
  mutate("AIC" = summary(MMR.MS.m3)$AICtab[1])

MMRDF = rbind(MMRDF1, MMRDF2, MMRDF3)

MMRDF$Response = "log Maximum Metabolic Rate" 
MMRDF$Response = "log Maximum Metabolic Rate" 
MMRDF$Response = "log Maximum Metabolic Rate" 

AASDF1 = subset(AAS.MS.Models, model == "Model 1") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 1")%>% 
  mutate("Marginal R2" = R2.Dataframe[7,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[7,2]) %>% 
  mutate("AIC" = summary(AAS.MS.m1)$AICtab[1])

AASDF2 = subset(AAS.MS.Models, model == "Model 2") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 2")%>% 
  mutate("Marginal R2" = R2.Dataframe[8,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[8,2]) %>% 
  mutate("AIC" = summary(AAS.MS.m2)$AICtab[1])

AASDF3 = subset(AAS.MS.Models, model == "Model 3") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 3")%>% 
  mutate("Marginal R2" = R2.Dataframe[9,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[9,2]) %>% 
  mutate("AIC" = summary(AAS.MS.m3)$AICtab[1])

AASDF = rbind(AASDF1, AASDF2, AASDF3)

AASDF$Response = "log Absolute Aerobic Scope" 
AASDF$Response = "log Absolute Aerobic Scope" 
AASDF$Response = "log Absolute Aerobic Scope" 

CTMDF = subset(CTM.Models, model == "Model 1") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 1")%>% 
  mutate("Marginal R2" = R2.Dataframe[10,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[10,2])%>%
  mutate("AIC" = summary(CTM.m1)$AICtab[1])

CTMDF$Response = "Critical Thermal Maximum" 

AllData = rbind(RMRDF, MMRDF, AASDF, CTMDF) 

AllData = AllData %>% select("Response", "Model Number", "term", "estimate", "std.error", "statistic", "df", 
                             "p.values", "AIC", 
                             "Marginal R2", "Conditional R2")
colnames(AllData) = c("Response Variable", "Model Number", "Model Parameter", "Estimate (beta)", "Std.Error", 
                      "T Statistic", "DF", "p", 
                      "AIC", "Marginal R2", "Conditional R2")

write.csv(AllData, "Model Output Results.csv")

###Testing Model Assumptions for Resting Metabic Rate: Model 1.

RMR.m1 = lmer(log(RMR) ~ log(Mass) + Mother.Acc + (1|Mother.ID) + (1|Father.ID), 
        REML = FALSE, data = RMR.NMS.Data)

summary(RMR.m1)$varcor 
##  Groups    Name        Std.Dev.  
##  Mother.ID (Intercept) 4.5652e-05
##  Father.ID (Intercept) 1.1354e-01
##  Residual              4.4331e-01
# Minimal variance explained by random effects, but retained for comparability.

    Res.Alone = ggplot(data = data.frame("Res" = residuals(RMR.m1, type = "deviance")), 
    aes(x = 1:nrow(RMR.NMS.Data), y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
    ylab("Deviance Residuals")

    Res.by.Fit = ggplot(data = data.frame("Res" = residuals(RMR.m1), "Fit" = predict(RMR.m1)), 
    aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
    ylab("Deviance Residuals")

    Res.by.Resp = ggplot(data = data.frame("Res" = residuals(RMR.m1), "Resp" = 
    log(RMR.NMS.Data$RMR)), aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
    xlab("Resting Metabolic Rate") + ylab("Deviance Residuals")

    QQNorm = ggplot(data.frame("Res" = residuals(RMR.m1), "Fit" = predict(RMR.m1)),
    aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = 
    "Resting Metabolic Rate Residuals")  

## Residuals likely non-normal. Slight high-end taper.

    qqPlot(residuals(RMR.m1, type = "deviance"), ylab = "Deviance Residuals")

## 83 81 
## 68 66
    # Within confidence intervals, however. Histogram to confirm, below.

    Resid.hist = ggplot(data.frame(Res = c(residuals(RMR.m1, type = "deviance")), Model = "Model1"), 
    aes(Res, fill = Model, colour = Model)) + geom_histogram(alpha = 0.5, bins = 20) + theme_bw() + 
    my.theme + xlab("Deviance Residuals") + ylab("Count") + scale_fill_manual(values = 
    c("cornflowerblue")) + scale_colour_manual(values = c("black")) + theme(legend.position = "none") + 
    geom_density(alpha = 0.5, aes(y=1.0/9 * ..count..))

    print(Resid.hist)

    # Reasonable distribution. Plotting residuals by predictors.

    Res.by.Mass = ggplot(data = data.frame("Mass" = RMR.NMS.Data$Mass,
    "Res" = residuals(RMR.m1)), aes(x = Mass, y = Res, colour = as.factor(RMR.NMS.Data$Mother.Acc))) + 
    theme_bw() + geom_point(size = 3, alpha = 0.5) + my.theme + 
    xlab("Mass (g)") + ylab("Deviance Residuals") + 
    scale_colour_manual(values = c("royalblue", "gold2")) + 
    guides(colour=guide_legend(title="Dam Acclimation \nTemperature"))

    Res.by.Dam = ggplot(data = data.frame("Dam" = as.factor(RMR.NMS.Data$Mother.Acc),
    "Res" = residuals(RMR.m1)), aes(x = Dam, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "royalblue") + my.theme + 
    xlab("Dam Acclimation Temperature (Degrees C)") + ylab("Deviance Residuals")
    
    grid.arrange(Res.by.Mass, Res.by.Dam, nrow = 2, top = 
    "Resting Metabolic Rate Residuals")  

    # Homoskedastic. Calculating p-values. 

    dwdata.RMR.NMS.m1 = data.frame(term = c(rownames(summary(RMR.m1)$coefficients)),
                       estimate = c(summary(RMR.m1)$coefficients[1:3]),
                       std.error = c(summary(RMR.m1)$coefficients[4:6]),
                       statistic = c(summary(RMR.m1)$coefficients[7:9]))

    p.values.m1 = c()
    for (i in 1:nrow(dwdata.RMR.NMS.m1)){
      p.values.m1[i] = pt(abs(dwdata.RMR.NMS.m1$statistic[i]), df = 8, lower.tail = FALSE)
    } 
    dwdata.RMR.NMS.m1$p.values = round(p.values.m1, 3)

    dwdata.RMR.NMS.m1$term = c("(Intercept)", "log Mass", "Dam Acclimation Temperature")

###Testing Model Assumptions for Resting Metabic Rate: Model 2.

RMR.m2 = lmer(log(RMR) ~ log(Mass) + Offspring.Acc + Mother.Acc + (1|Mother.ID) + (1|Father.ID), 
        REML = FALSE, data = RMR.NMS.Data)

summary(RMR.m2)$varcor 
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.00000 
##  Father.ID (Intercept) 0.12491 
##  Residual              0.44152
# Maternal ID alone explaining zero variance.

    Res.Alone = ggplot(data = data.frame("Res" = residuals(RMR.m2, type = "deviance")), 
    aes(x = 1:nrow(RMR.NMS.Data), y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
    ylab("Deviance Residuals")

    Res.by.Fit = ggplot(data = data.frame("Res" = residuals(RMR.m2, type = "deviance"), 
    "Fit" = predict(RMR.m2)), aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
    ylab("Deviance Residuals")

    Res.by.Resp = ggplot(data = data.frame("Res" = residuals(RMR.m2, type = "deviance"), 
    "Resp" = log(RMR.NMS.Data$RMR)), aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
    xlab("Resting Metabolic Rate") + ylab("Deviance Residuals")

    QQNorm = ggplot(data.frame("Res" = residuals(RMR.m2, type = "deviance"), 
    "Fit" = predict(RMR.m2)), aes(sample=Res)) + stat_qq(colour = "gold") + 
    stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = 
    "Resting Metabolic Rate Residuals")  

## Residuals similar to model 1. 

    qqPlot(residuals(RMR.m2, type = "deviance"), ylab = "Deviance Residuals")

## 83 81 
## 68 66
    # Again, generally within confidence intervals.

    Res.by.Mass = ggplot(data = data.frame("Mass" = RMR.NMS.Data$Mass,
    "Res" = residuals(RMR.m2)), aes(x = Mass, y = Res, colour = as.factor(RMR.NMS.Data$Mother.Acc))) + 
    theme_bw() + geom_point(size = 3, alpha = 0.5) + my.theme + 
    xlab("Mass (g)") + ylab("Deviance Residuals") + 
    scale_colour_manual(values = c("royalblue", "gold2")) + 
    guides(colour=guide_legend(title="Dam Acclimation \nTemperature"))

    Res.by.Dam = ggplot(data = data.frame("Dam" = as.factor(RMR.NMS.Data$Mother.Acc),
    "Res" = residuals(RMR.m2)), aes(x = Dam, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "royalblue") + my.theme.diag + 
    xlab("Dam Acclimation Temperature \n (Degrees C)") + ylab("Deviance Residuals")
    
    Res.by.Off = ggplot(data = data.frame("Offspring" = as.factor(RMR.NMS.Data$Offspring.Acc),
    "Res" = residuals(RMR.m2)), aes(x = Offspring, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "gold2") + my.theme.diag +
    xlab("Dam Acclimation Temperature \n (Degrees C)") + ylab("Deviance Residuals")

    layout <- rbind(c(1,1),
             c(2,3))

    grid.arrange(Res.by.Mass, Res.by.Dam, Res.by.Off, nrow = 2, top = 
    "Resting Metabolic Rate Residuals", layout_matrix = layout)  

    # Again, homoskedastic. Calculating p-values.

    dwdata.RMR.NMS.m2 = data.frame(term = c(rownames(summary(RMR.m2)$coefficients)),
                       estimate = c(summary(RMR.m2)$coefficients[1:4]),
                       std.error = c(summary(RMR.m2)$coefficients[5:8]),
                       statistic = c(summary(RMR.m2)$coefficients[9:12]))

    p.values.m2 = c()
    for (i in 1:nrow(dwdata.RMR.NMS.m2)){
      p.values.m2[i] = pt(abs(dwdata.RMR.NMS.m2$statistic[i]), df = 8, lower.tail = FALSE)
    } 
    dwdata.RMR.NMS.m2$p.values = round(p.values.m2, 3)

    dwdata.RMR.NMS.m2$term = c("(Intercept)", "log Mass", "Offspring Acclimation Temperature", 
    "Dam Acclimation Temperature")

###Testing Model Assumptions for Resting Metabic Rate: Model 3.

RMR.m3 = lmer(log(RMR) ~ log(Mass) + Offspring.Acc + Mother.Acc + Offspring.Acc:Mother.Acc + 
        (1|Mother.ID) + (1|Father.ID), REML = FALSE, data = RMR.NMS.Data)

summary(RMR.m3)$varcor 
##  Groups    Name        Std.Dev. 
##  Mother.ID (Intercept) 0.0083944
##  Father.ID (Intercept) 0.1284146
##  Residual              0.4390322
# Maternal ID alone explaining zero variance.

    Res.Alone = ggplot(data = data.frame("Res" = residuals(RMR.m3, type = "deviance")), 
    aes(x = 1:nrow(RMR.NMS.Data), y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
    ylab("Deviance Residuals")

    Res.by.Fit = ggplot(data = data.frame("Res" = residuals(RMR.m3, type = "deviance"), 
    "Fit" = predict(RMR.m3)), aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
    ylab("Deviance Residuals")

    Res.by.Resp = ggplot(data = data.frame("Res" = residuals(RMR.m3, type = "deviance"), 
    "Resp" = log(RMR.NMS.Data$RMR)), aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
    xlab("Resting Metabolic Rate") + ylab("Deviance Residuals")

    QQNorm = ggplot(data.frame("Res" = residuals(RMR.m3, type = "deviance"), 
    "Fit" = predict(RMR.m3)), aes(sample=Res)) + stat_qq(colour = "gold") + 
    stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = 
    "Resting Metabolic Rate Residuals")  

## Again, slight tails. 

    qqPlot(residuals(RMR.m3, type = "deviance"), ylab = "Deviance Residuals")

## 83 81 
## 68 66
## Largely within confidence intervals.

    Resid.hist = ggplot(data.frame(Res = c(residuals(RMR.m3, type = "deviance")), Model = "Model1"), 
    aes(Res, fill = Model, colour = Model)) + geom_histogram(alpha = 0.5, bins = 20) + theme_bw() + 
    my.theme + xlab("Deviance Residuals") + ylab("Count") + scale_fill_manual(values = 
    c("firebrick")) + scale_colour_manual(values = c("black")) + theme(legend.position = "none") + 
    geom_density(alpha = 0.5, aes(y=0.5/5 * ..count..))

    print(Resid.hist)

    # Slightly bimodal. Assessing predictors.

    Res.by.Mass = ggplot(data = data.frame("Mass" = RMR.NMS.Data$Mass,
    "Res" = residuals(RMR.m2)), aes(x = Mass, y = Res, colour = as.factor(RMR.NMS.Data$Mother.Acc))) + 
    theme_bw() + geom_point(size = 3, alpha = 0.5) + my.theme.diag + 
    xlab("Mass (g)") + ylab("Deviance Residuals") + 
    scale_colour_manual(values = c("royalblue", "gold2")) + 
    guides(colour=guide_legend(title="Dam Acclimation Temperature")) + 
    theme(legend.position = "top")

    Res.by.Dam = ggplot(data = data.frame("Dam" = as.factor(RMR.NMS.Data$Mother.Acc),
    "Res" = residuals(RMR.m2)), aes(x = Dam, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "royalblue") + my.theme.diag + 
    xlab("Dam Acclimation Temperature \n (Degrees C)") + ylab("Deviance Residuals")
    
    Res.by.Off = ggplot(data = data.frame("Offspring" = as.factor(RMR.NMS.Data$Offspring.Acc),
    "Res" = residuals(RMR.m2)), aes(x = Offspring, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "gold2") + my.theme.diag + 
    xlab("Dam Acclimation Temperature \n (Degrees C)") + ylab("Deviance Residuals")

    Off.by.Dam = ggplot(data = data.frame("Dam" = as.factor(RMR.NMS.Data$Mother.Acc),
    "Res" = residuals(RMR.m2)), aes(x = Dam, y = Res, fill = as.factor(RMR.NMS.Data$Offspring.Acc))) + 
    theme_bw() + geom_boxplot() + my.theme.diag + scale_fill_manual(values = c("firebrick", "coral")) +
    xlab("Dam Acclimation Temperature (Degrees C)") + ylab("Deviance Residuals") + 
    guides(fill=guide_legend(title="Offspring Acclimation \n Temperature")) + 
    theme(legend.position = "top")

    layout <- rbind(c(1,1,2,3),
             c(1,1,4,4))

    grid.arrange(Res.by.Mass, Res.by.Dam, Res.by.Off, Off.by.Dam, nrow = 2, top = 
    "Resting Metabolic Rate Residuals", layout_matrix = layout)  

    # Largely homoskedastic. Calculating p-values.

    dwdata.RMR.NMS.m3 = data.frame(term = c(rownames(summary(RMR.m3)$coefficients)),
                       estimate = c(summary(RMR.m3)$coefficients[1:5]),
                       std.error = c(summary(RMR.m3)$coefficients[6:10]),
                       statistic = c(summary(RMR.m3)$coefficients[11:15]))

    p.values.m3 = c()
    for (i in 1:nrow(dwdata.RMR.NMS.m3)){
      p.values.m3[i] = pt(abs(dwdata.RMR.NMS.m3$statistic[i]), df = 8, lower.tail = FALSE)
    } 
    dwdata.RMR.NMS.m3$p.values = round(p.values.m3, 3)

    dwdata.RMR.NMS.m3$term = c("(Intercept)", "log Mass", "Offspring Acclimation Temperature", 
    "Dam Acclimation Temperature", "Dam Acclimation Temperature :\n Offspring Acclimation Temperature" )

## Comparing models by plot.

dwdata.RMR.NMS.m1 = dwdata.RMR.NMS.m1 %>% mutate(model = "Model 1")
dwdata.RMR.NMS.m2 = dwdata.RMR.NMS.m2 %>% mutate(model = "Model 2")
dwdata.RMR.NMS.m3 = dwdata.RMR.NMS.m3 %>% mutate(model = "Model 3")

RMR.NMS.Models = rbind(dwdata.RMR.NMS.m1, dwdata.RMR.NMS.m2, dwdata.RMR.NMS.m3)

RMR.NMS.Models$estimate = RMR.NMS.Models$estimate*5
RMR.NMS.Models$std.error = RMR.NMS.Models$std.error*5 

RMR.NMS.Models %>%
  dwplot(style = "distribution", show_intercept = FALSE) + theme_bw() + 
  theme(legend.justification = c(-0.1,-0.1), legend.position = c(-0.8, 0.01)) + my.theme.dw + 
  xlab("Coefficient") + ylab("") + geom_vline(xintercept = 0, colour = "black", linetype = 2) + 
  scale_fill_manual(values = c("mediumseagreen", "mediumorchid", "cornflowerblue")) + 
  scale_colour_manual(values = c("black", "black", "black"))  + scale_x_continuous(limits = c(-5.0,5.0), 
  breaks = c(-5.0, -2.5, 0, 2.5, 5.0), label = c("-1.00", "-0.50", "0", "0.50", "1.00"))

ggsave("/home/joshk/git_repositories/BrookTrout_TG/Figures/log Resting Metabolic Rate Model Comparisons.jpeg", dpi = 800, limitsize = TRUE)

RMR.NMS.Models$estimate = RMR.NMS.Models$estimate/5
RMR.NMS.Models$std.error = RMR.NMS.Models$std.error/5 

RMR.NMS.Fit.m1 = data.frame("Response" = c(log(RMR.NMS.Data$RMR)), "Fit" = c(predict(RMR.m1)))
RMR.NMS.Fit.m2 = data.frame("Response" = c(log(RMR.NMS.Data$RMR)), "Fit" = c(predict(RMR.m2)))
RMR.NMS.Fit.m3 = data.frame("Response" = c(log(RMR.NMS.Data$RMR)), "Fit" = c(predict(RMR.m3)))
RMR.NMS.Fit.m1 = RMR.NMS.Fit.m1 %>% mutate(model = "Model 1")
RMR.NMS.Fit.m2 = RMR.NMS.Fit.m2 %>% mutate(model = "Model 2")
RMR.NMS.Fit.m3 = RMR.NMS.Fit.m3 %>% mutate(model = "Model 3")

RMR.NMS.Fit.All = rbind(RMR.NMS.Fit.m1, RMR.NMS.Fit.m2, RMR.NMS.Fit.m3)
RMR.NMS.Fits = ggplot(data = RMR.NMS.Fit.All, aes(x = Fit, y = Response, fill = model, colour = model,
  linetype = model)) + theme_bw() + my.theme + geom_point(pch=21, size = 3, alpha = 0.5) + 
  geom_smooth(method = "lm") + xlab("Predicted log Resting Metabolic Rate \n (log mg O2/h)") + 
  ylab("Measured log Resting \n Metabolic Rate \n (log mg O2/h)") + 
  scale_fill_manual(values = c("mediumseagreen", 
  "mediumorchid", "cornflowerblue")) + scale_colour_manual(values = c("black", 
  "black", "black")) + theme(legend.title=element_blank())

print(RMR.NMS.Fits)

ggsave("/home/joshk/git_repositories/BrookTrout_TG/Figures/log Resting Metabolic Rate Model Fit Comparisons.jpeg", dpi = 800, limitsize = TRUE)

###Testing Model Assumptions for Maximum Metabic Rate: Model 1.

    MMR.m1 = lmer(log(MMR) ~ log(Mass) + (1|Mother.ID) + (1|Father.ID), 
        REML = FALSE, data = MMR.NMS.Data)

    summary(MMR.m1)$varcor
##  Groups    Name        Std.Dev.  
##  Mother.ID (Intercept) 0.00054796
##  Father.ID (Intercept) 0.01631463
##  Residual              0.25085402
    # Maternal ID explains zero residual variance. 

    Res.Alone = ggplot(data = data.frame("Res" = residuals(MMR.m1, type = "deviance")), 
    aes(x = 1:nrow(MMR.NMS.Data), y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
    ylab("Deviance Residuals")

    Res.by.Fit = ggplot(data = data.frame("Res" = residuals(MMR.m1, type = "deviance"), 
    "Fit" = predict(MMR.m1)), aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
    ylab("Deviance Residuals")

    Res.by.Resp = ggplot(data = data.frame("Res" = residuals(MMR.m1, type = "deviance"), 
    "Resp" = MMR.NMS.Data$MMR), aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
    xlab("Resting Metabolic Rate") + ylab("Deviance Residuals")

    QQNorm = ggplot(data.frame("Res" = residuals(MMR.m1, type = "deviance"), 
    "Fit" = predict(MMR.m1)), aes(sample=Res)) + stat_qq(colour = "gold") + 
    stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = 
    "Maximum Metabolic Rate Residuals")  

    # Wow! One very destinct outlier. Assessing individual.

    MMR.NMS.Data[which(residuals(MMR.m1, type = "deviance") < -2.0),]
##    Offspring.Acc Mother.Acc Father.Acc Mother.ID Father.ID Mass    RMR    MMR
## 83             1          1          1         3         4    3 0.1118 0.1165
##        AAS  CTM Parent.group Tank Tank.Box Family
## 83 0.05248 28.3            1    3       3C     15
    # Remarkably low maximum metabolic rate and absolute aerobic scope. Removing from both datasets.

    MMR.NMS.Data = MMR.NMS.Data[-c(which(residuals(MMR.m1, type = "deviance") < -2.0)),]
    AAS.NMS.Data = AAS.NMS.Data[-c(which(AAS.NMS.Data$AAS == 0.05248)),]

Recalculting AICc for Maximum Metabolic Rate - Model Iterations.

    MMR.model = lmer(log(MMR) ~ Mother.Acc*Father.Acc*Offspring.Acc + log(Mass) + (1|Mother.ID) +(1|Father.ID), 
    REML = FALSE, data = MMR.NMS.Data)  

    AIC.Results.MMR = subset(dredge(MMR.model, rank = "AICc", 
        extra = list(R2 = function(x) r.squaredGLMM(x))), delta < 5)

Select.Models.MMR = c()
Delta.MMR = c()
Response = c()
Marg.R2 = c() 
Cond.R2 = c()

  for (i in 1:length(attr(AIC.Results.MMR, "model.calls"))){
    Select.Models.MMR[i] = gsub(".*~ ", "", paste(attr(AIC.Results.MMR, "model.calls")[[i]])[2])
    Delta.MMR[i] = AIC.Results.MMR$delta[i]
    Response[i] = gsub(" ~.*", "", paste(attr(AIC.Results.MMR, "model.calls")[[i]])[2])
    Marg.R2[i] = AIC.Results.MMR$R21[i]
    Cond.R2[i] = AIC.Results.MMR$R22[i]
    if (i == length(attr(AIC.Results.MMR, "model.calls"))){
      Table = data.frame("Reponse.Type" = "Mass Independant", "Response" = Response, 
      "Delta.AIC" = Delta.MMR, "Marginal.R2" = Marg.R2, "Conditional.R2" = Cond.R2,
      "Equation" = Select.Models.MMR)
      assign(paste("AIC.Selection", "MMR", sep = "_"), Table)
    }
  }

    AIC.Selection.NMS = rbind(AIC.Selection_NMS_1, AIC.Selection_MMR, AIC.Selection_NMS_3, AIC.Selection_NMS_4)

    AIC.Selection = rbind(AIC.Selection.MS, AIC.Selection.NMS)
    print(AIC.Selection, right = F)
##    Reponse.Type     Response    Delta.AIC Marginal.R2 Conditional.R2
## 1  Mass-Specific    log(RMR.MS) 0.0000000 0.09689341  0.23258538    
## 2  Mass-Specific    log(RMR.MS) 0.4320103 0.05317331  0.24043771    
## 3  Mass-Specific    log(RMR.MS) 1.2633111 0.11949983  0.22942058    
## 4  Mass-Specific    log(RMR.MS) 1.6197476 0.07652765  0.23332101    
## 5  Mass-Specific    log(RMR.MS) 2.1592505 0.09688383  0.23257954    
## 6  Mass-Specific    log(RMR.MS) 3.4158810 0.11948859  0.22852419    
## 7  Mass-Specific    log(RMR.MS) 3.4339163 0.11966268  0.23078182    
## 8  Mass-Specific    log(RMR.MS) 3.4465742 0.11949854  0.22941555    
## 9  Mass-Specific    log(RMR.MS) 3.7625995 0.07663967  0.23491827    
## 10 Mass-Specific    log(MMR.MS) 0.0000000 0.12081940  0.14201563    
## 11 Mass-Specific    log(MMR.MS) 1.1337870 0.12482063  0.14707607    
## 12 Mass-Specific    log(MMR.MS) 1.6198084 0.07304761  0.14942805    
## 13 Mass-Specific    log(MMR.MS) 1.6632302 0.12283706  0.14581838    
## 14 Mass-Specific    log(MMR.MS) 2.7260488 0.07701051  0.15563521    
## 15 Mass-Specific    log(MMR.MS) 2.7981663 0.12640953  0.14714116    
## 16 Mass-Specific    log(MMR.MS) 2.8091275 0.12687437  0.15094298    
## 17 Mass-Specific    log(MMR.MS) 3.1451250 0.12529971  0.14712299    
## 18 Mass-Specific    log(MMR.MS) 4.4679396 0.07829618  0.15562601    
## 19 Mass-Specific    log(MMR.MS) 4.5445400 0.12823397  0.15073801    
## 20 Mass-Specific    log(MMR.MS) 4.7897724 0.12703622  0.14722626    
## 21 Mass-Specific    log(MMR.MS) 4.8318359 0.12737623  0.15100909    
## 22 Mass-Specific    log(AAS.MS) 0.0000000 0.06030224  0.06030224    
## 23 Mass-Specific    log(AAS.MS) 1.6058206 0.06331924  0.06331924    
## 24 Mass-Specific    log(AAS.MS) 1.6547206 0.06304827  0.06304827    
## 25 Mass-Specific    log(AAS.MS) 2.7150419 0.06921105  0.06921105    
## 26 Mass-Specific    log(AAS.MS) 3.2960644 0.06600676  0.06600676    
## 27 Mass-Specific    log(AAS.MS) 3.8171171 0.06312377  0.06312377    
## 28 Mass-Specific    log(AAS.MS) 4.4085037 0.07201229  0.07201229    
## 29 Mass Independant log(RMR)    0.0000000 0.13398490  0.18729777    
## 30 Mass Independant log(RMR)    1.2602254 0.14177631  0.20537453    
## 31 Mass Independant log(RMR)    1.6073708 0.15042976  0.21764843    
## 32 Mass Independant log(RMR)    1.6474074 0.12869909  0.17308051    
## 33 Mass Independant log(RMR)    2.7321064 0.13274864  0.17071098    
## 34 Mass Independant log(RMR)    2.8878776 0.13578794  0.19007496    
## 35 Mass Independant log(RMR)    3.2743628 0.14463523  0.20289714    
## 36 Mass Independant log(RMR)    4.0588194 0.14000948  0.18880081    
## 37 Mass Independant log(RMR)    4.0950899 0.09988873  0.21612058    
## 38 Mass Independant log(RMR)    4.5640317 0.14806374  0.20017371    
## 39 Mass Independant log(RMR)    4.8488012 0.13927880  0.19821267    
## 40 Mass Independant log(MMR)    0.0000000 0.51073885  0.52414430    
## 41 Mass Independant log(MMR)    1.1519770 0.50259058  0.51232882    
## 42 Mass Independant log(MMR)    1.6571485 0.51244418  0.52665944    
## 43 Mass Independant log(MMR)    2.3452891 0.51490234  0.52733792    
## 44 Mass Independant log(MMR)    2.7795841 0.50375882  0.51398621    
## 45 Mass Independant log(MMR)    3.2968819 0.50262005  0.51249178    
## 46 Mass Independant log(MMR)    3.4572491 0.50709834  0.51576702    
## 47 Mass Independant log(MMR)    4.6305941 0.48890967  0.49534628    
## 48 Mass Independant log(MMR)    4.8358232 0.50402584  0.51443503    
## 49 Mass Independant log(MMR)    4.9419345 0.50379582  0.51418385    
## 50 Mass Independant AAS         0.0000000 0.37930675  0.37930675    
## 51 Mass Independant AAS         1.5094369 0.38174466  0.38174466    
## 52 Mass Independant AAS         2.1236991 0.37949846  0.37949846    
## 53 Mass Independant AAS         2.6531066 0.38559940  0.38559940    
## 54 Mass Independant AAS         3.6519021 0.38196577  0.38196577    
## 55 Mass Independant AAS         4.3197068 0.37952426  0.37952426    
## 56 Mass Independant AAS         4.8437724 0.38574420  0.38574420    
## 57 Mass Independant CTM         0.0000000 0.46835156  0.50879302    
##    Equation                                                                                                                                  
## 1  Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 2  Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                         
## 3  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                               
## 4  Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 5  Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                 
## 6  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                       
## 7  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                    
## 8  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                    
## 9  Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                 
## 10 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 11 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                               
## 12 Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                         
## 13 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                 
## 14 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 15 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                    
## 16 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                    
## 17 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                       
## 18 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                 
## 19 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc + Mother.Acc:Offspring.Acc         
## 20 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc + Mother.Acc:Offspring.Acc            
## 21 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc + Father.Acc:Offspring.Acc            
## 22 Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                         
## 23 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 24 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 25 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                 
## 26 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                               
## 27 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                 
## 28 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                    
## 29 log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                
## 30 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 31 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                     
## 32 Father.Acc + log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                   
## 33 Father.Acc + log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                           
## 34 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                   
## 35 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                        
## 36 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                           
## 37 log(Mass) + (1 | Mother.ID) + (1 | Father.ID)                                                                                             
## 38 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc + Mother.Acc:Offspring.Acc
## 39 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                        
## 40 log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                             
## 41 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 42 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 43 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                     
## 44 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                   
## 45 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                     
## 46 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                        
## 47 log(Mass) + (1 | Mother.ID) + (1 | Father.ID)                                                                                             
## 48 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                           
## 49 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                        
## 50 log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                             
## 51 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 52 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 53 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                     
## 54 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                   
## 55 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                     
## 56 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                        
## 57 Offspring.Acc
    write.csv(AIC.Selection, "Model_Outputs_AICc.csv", row.names = FALSE)

###Reassessing Model Assumptions for Maximum Metabic Rate: Model 1.

    MMR.m1 = lmer(log(MMR) ~ log(Mass) + Offspring.Acc + (1|Mother.ID) + (1|Father.ID), 
        REML = FALSE, data = MMR.NMS.Data)

    summary(MMR.m1)$AICtab[1]
##      AIC 
## -82.4201
    # AIC dramatically improved.

    summary(MMR.m1)$varcor
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.000000
##  Father.ID (Intercept) 0.032587
##  Residual              0.194153
    # Maternal ID again explaining no residual variance.

    Res.Alone = ggplot(data = data.frame("Res" = residuals(MMR.m1, type = "deviance")), 
    aes(x = 1:nrow(MMR.NMS.Data), y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
    ylab("Deviance Residuals")

    Res.by.Fit = ggplot(data = data.frame("Res" = residuals(MMR.m1, type = "deviance"), 
    "Fit" = predict(MMR.m1)), aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
    ylab("Deviance Residuals")

    Res.by.Resp = ggplot(data = data.frame("Res" = residuals(MMR.m1, type = "deviance"), 
    "Resp" = MMR.NMS.Data$MMR), aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
    xlab("Resting Metabolic Rate") + ylab("Deviance Residuals")

    QQNorm = ggplot(data.frame("Res" = residuals(MMR.m1, type = "deviance"), 
    "Fit" = predict(MMR.m1)), aes(sample=Res)) + stat_qq(colour = "gold") + 
    stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = 
    "Maximum Metabolic Rate Residuals")  

    qqPlot(residuals(MMR.m1, type = "deviance"), ylab = "Deviance Residuals")

## 198  47 
## 194  47
    # Residuals horribly non-normal.

    Resid.hist = ggplot(data.frame(Res = c(residuals(MMR.m1, type = "deviance")), Model = "Model1"), 
    aes(Res, fill = Model, colour = Model)) + geom_histogram(alpha = 0.5, bins = 20) + theme_bw() + 
    my.theme + xlab("Deviance Residuals") + ylab("Count") + scale_fill_manual(values = 
    c("black")) + scale_colour_manual(values = c("navajowhite")) + theme(legend.position = "none") + 
    geom_density(alpha = 0.5, aes(y=0.2/5 * ..count..))
    
    print(Resid.hist)

    # Bimodal. Assessing predictors by residuals.

    Res.by.Mass = ggplot(data = data.frame("Mass" = RMR.NMS.Data$Mass,
    "Res" = residuals(RMR.m2)), aes(x = Mass, y = Res, colour = as.factor(RMR.NMS.Data$Offspring.Acc))) + 
    theme_bw() + geom_point(size = 3, alpha = 0.5) + my.theme.diag + 
    xlab("Mass (g)") + ylab("Deviance Residuals") + 
    scale_colour_manual(values = c("royalblue", "gold2")) + 
    guides(colour=guide_legend(title="Offspring Acclimation Temperature")) + 
    theme(legend.position = "top")

    Res.by.Off = ggplot(data = data.frame("Off" = as.factor(RMR.NMS.Data$Offspring.Acc),
    "Res" = residuals(RMR.m2)), aes(x = Off, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "royalblue") + my.theme.diag + 
    xlab("Offspring Acclimation Temperature (Degrees C)") + ylab("Deviance Residuals")

    QQNorm.Diag = ggplot(data.frame("Res" = residuals(MMR.m1, type = "deviance"), 
    "Fit" = predict(MMR.m1), "Off" = as.factor(MMR.NMS.Data$Offspring.Acc)), 
    aes(sample=Res, colour = Off)) + stat_qq() + 
    stat_qq_line() + my.theme + theme_bw() + 
    scale_colour_manual(values = c("black", "thistle")) + theme(legend.position = "top") +
    guides(colour=guide_legend(title="Offspring Acclimation Temperature"))

    layout <- rbind(c(1,1,2,2),
             c(NA,3,3,NA))

    grid.arrange(Res.by.Mass, Res.by.Off, QQNorm.Diag, nrow = 2, top = 
    "Maximum Metabolic Rate Residuals", layout_matrix = layout)  

    # No clear heteroskedasticity, nor difference in sample size. Assessing a histogram.

    Resid.hist = ggplot(data.frame(Res = c(residuals(MMR.m1, type = "deviance")), 
    "Offspring" = as.factor(MMR.NMS.Data$Offspring.Acc)), 
    aes(Res, fill = Offspring, colour = Offspring)) + geom_histogram(alpha = 0.5, bins = 20) + theme_bw() + 
    my.theme + xlab("Deviance Residuals") + ylab("Count") + scale_fill_manual(values = 
    c("black", "mediumorchid")) + scale_colour_manual(values = c("white", "black")) + 
    theme(legend.position = "none") + 
    geom_density(alpha = 0.5, aes(y=0.2/4 * ..count..))
    
    print(Resid.hist)

    ## Equivalent data distributions. Continuing on.

    dwdata.MMR.NMS.m1 = data.frame(term = c(rownames(summary(MMR.m1)$coefficients)),
                       estimate = c(summary(MMR.m1)$coefficients[1:3]),
                       std.error = c(summary(MMR.m1)$coefficients[4:6]),
                       statistic = c(summary(MMR.m1)$coefficients[7:9]))

    p.values.m1 = c()
    for (i in 1:nrow(dwdata.MMR.NMS.m1)){
      p.values.m1[i] = pt(abs(dwdata.MMR.NMS.m1$statistic[i]), df = 8, lower.tail = FALSE)
    } 
    dwdata.MMR.NMS.m1$p.values = round(p.values.m1, 3)

    dwdata.MMR.NMS.m1$term = c("(Intercept)", "log Mass", "Offspring Acclimation Temperature")

###Testing Model Assumptions for Maximum Metabic Rate: Model 2.

    MMR.m2 = lmer(log(MMR) ~ log(Mass) + Offspring.Acc + Father.Acc +  
        (1|Mother.ID) + (1|Father.ID), REML = FALSE, data = MMR.NMS.Data)

    summary(MMR.m2)$varcor
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.000000
##  Father.ID (Intercept) 0.027424
##  Residual              0.194068
    # Maternal and paternal ID explain very little residual variance. 

    Res.Alone = ggplot(data = data.frame("Res" = residuals(MMR.m2, type = "deviance")), 
    aes(x = 1:nrow(MMR.NMS.Data), y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
    ylab("Deviance Residuals")

    Res.by.Fit = ggplot(data = data.frame("Res" = residuals(MMR.m2, type = "deviance"), 
    "Fit" = predict(MMR.m2)), aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
    ylab("Deviance Residuals")

    Res.by.Resp = ggplot(data = data.frame("Res" = residuals(MMR.m2, type = "deviance"), 
    "Resp" = MMR.NMS.Data$MMR), aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
    xlab("Resting Metabolic Rate") + ylab("Deviance Residuals")

    QQNorm = ggplot(data.frame("Res" = residuals(MMR.m2, type = "deviance"), 
    "Fit" = predict(MMR.m2)), aes(sample=Res)) + stat_qq(colour = "gold") + 
    stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = 
    "Maximum Metabolic Rate Residuals")  

    # Normality far from met again.

    qqPlot(residuals(MMR.m2, type = "deviance"), ylab = "Deviance Residuals")

## 198  47 
## 194  47
    Resid.hist = ggplot(data.frame(Res = c(residuals(MMR.m2, type = "deviance")), Model = "Model1"), 
    aes(Res, fill = Model, colour = Model)) + geom_histogram(alpha = 0.5, bins = 20) + theme_bw() + 
    my.theme + xlab("Deviance Residuals") + ylab("Count") + scale_fill_manual(values = 
    c("cornflowerblue")) + scale_colour_manual(values = c("black")) + theme(legend.position = "none") + 
    geom_density(alpha = 0.5, aes(y=0.2/5 * ..count..))
    
    print(Resid.hist)

    # Arguably more pronounced bimodality. Assessing parameter-specific variance.

    Res.by.Mass = ggplot(data = data.frame("Mass" = MMR.NMS.Data$Mass,
    "Res" = residuals(MMR.m2, type = "deviance")), aes(x = Mass, y = Res, 
    colour = as.factor(MMR.NMS.Data$Offspring.Acc))) + 
    theme_bw() + geom_point(size = 3, alpha = 0.5) + my.theme.diag + 
    xlab("Mass (g)") + ylab("Deviance Residuals") + 
    scale_colour_manual(values = c("royalblue", "gold2")) + 
    guides(colour=guide_legend(title="Offspring Acclimation Temperature")) + 
    theme(legend.position = "top")

    Res.by.Off = ggplot(data = data.frame("Off" = as.factor(MMR.NMS.Data$Offspring.Acc),
    "Res" = residuals(MMR.m2, type = "deviance")), aes(x = Off, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "royalblue") + my.theme.diag + 
    xlab("Offspring Acclimation Temperature (Degrees C)") + ylab("Deviance Residuals")

    Res.by.Pat = ggplot(data = data.frame("Pat" = as.factor(MMR.NMS.Data$Father.Acc),
    "Res" = residuals(MMR.m2, type = "deviance")), aes(x = Pat, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "mediumseagreen") + my.theme.diag + 
    xlab("Sire Acclimation Temperature (Degrees C)") + ylab("Deviance Residuals")

    QQNorm.by.Off = ggplot(data.frame("Res" = residuals(MMR.m2, type = "deviance"), 
    "Fit" = predict(MMR.m2), "Off" = as.factor(MMR.NMS.Data$Offspring.Acc)), 
    aes(sample=Res, colour = Off)) + stat_qq() + 
    stat_qq_line() + my.theme + theme_bw() + 
    scale_colour_manual(values = c("mediumseagreen", "mediumorchid")) + theme(legend.position = "top") +
    guides(colour=guide_legend(title="Offspring Acclimation Temperature"))

    QQNorm.by.Pat = ggplot(data.frame("Res" = residuals(MMR.m2, type = "deviance"), 
    "Fit" = predict(MMR.m2), "Pat" = as.factor(MMR.NMS.Data$Father.Acc)), 
    aes(sample=Res, colour = Pat)) + stat_qq() + 
    stat_qq_line() + my.theme + theme_bw() + 
    scale_colour_manual(values = c("royalblue", "gold2")) + theme(legend.position = "top") +
    guides(colour=guide_legend(title="Sire Acclimation Temperature"))

    grid.arrange(Res.by.Mass, Res.by.Off, Res.by.Pat, QQNorm.by.Off, QQNorm.by.Pat,
    nrow = 3, top = "Maximum Metabolic Rate Residuals") 

    # Perhaps a slight heterogeneity across offspring acclimation temperature, but very minimal.

    # Calculating p-values.

    dwdata.MMR.NMS.m2 = data.frame(term = c(rownames(summary(MMR.m2)$coefficients)),
                       estimate = c(summary(MMR.m2)$coefficients[1:4]),
                       std.error = c(summary(MMR.m2)$coefficients[5:8]),
                       statistic = c(summary(MMR.m2)$coefficients[9:12]))

    p.values.m2 = c()
    for (i in 1:nrow(dwdata.MMR.NMS.m2)){
      p.values.m2[i] = pt(abs(dwdata.MMR.NMS.m2$statistic[i]), df = 8, lower.tail = FALSE)
    } 
    dwdata.MMR.NMS.m2$p.values = round(p.values.m2, 3)

    dwdata.MMR.NMS.m2$term = c("(Intercept)", "log Mass", "Offspring Acclimation Temperature",
        "Sire Acclimation Temperature")

###Testing Model Assumptions for Maximum Metabic Rate: Model 3.

    MMR.m3 = lmer(log(MMR) ~ log(Mass) + Offspring.Acc + Mother.Acc + 
    (1|Mother.ID) + (1|Father.ID), REML = FALSE, data = MMR.NMS.Data)

    summary(MMR.m3)$varcor
##  Groups    Name        Std.Dev.  
##  Mother.ID (Intercept) 9.5575e-10
##  Father.ID (Intercept) 3.3598e-02
##  Residual              1.9388e-01
    # Paternal ID alone explaining residual variance. 

    Res.Alone = ggplot(data = data.frame("Res" = residuals(MMR.m3, type = "deviance")), 
    aes(x = 1:nrow(MMR.NMS.Data), y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
    ylab("Deviance Residuals")

    Res.by.Fit = ggplot(data = data.frame("Res" = residuals(MMR.m3), "Fit" = predict(MMR.m3)), 
    aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
    ylab("Deviance Residuals")

    Res.by.Resp = ggplot(data = data.frame("Res" = residuals(MMR.m3), "Resp" = 
    MMR.NMS.Data$MMR), aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
    xlab("Resting Metabolic Rate") + ylab("Deviance Residuals")

    QQNorm = ggplot(data.frame("Res" = residuals(MMR.m3), "Fit" = predict(MMR.m2)),
    aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = 
    "Maximum Metabolic Rate Residuals")  

    # Similar patterns to previous two models.

    Resid.hist = ggplot(data.frame(Res = c(residuals(MMR.m3, type = "deviance")), Model = "Model1"), 
    aes(Res, fill = Model, colour = Model)) + geom_histogram(alpha = 0.5, bins = 20) + theme_bw() + 
    my.theme + xlab("Deviance Residuals") + ylab("Count") + scale_fill_manual(values = 
    c("cornflowerblue")) + scale_colour_manual(values = c("black")) + theme(legend.position = "none") + 
    geom_density(alpha = 0.5, aes(y=0.2/5 * ..count..))
    
    print(Resid.hist)

    Res.by.Mass = ggplot(data = data.frame("Mass" = MMR.NMS.Data$Mass,
    "Res" = residuals(MMR.m3, type = "deviance")), aes(x = Mass, y = Res, 
    colour = as.factor(MMR.NMS.Data$Offspring.Acc))) + 
    theme_bw() + geom_point(size = 3, alpha = 0.5) + my.theme.diag + 
    xlab("Mass (g)") + ylab("Deviance Residuals") + 
    scale_colour_manual(values = c("royalblue", "gold2")) + 
    guides(colour=guide_legend(title="Offspring Acclimation Temperature")) + 
    theme(legend.position = "top")

    Res.by.Off = ggplot(data = data.frame("Off" = as.factor(MMR.NMS.Data$Offspring.Acc),
    "Res" = residuals(MMR.m3, type = "deviance")), aes(x = Off, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "royalblue") + my.theme.diag + 
    xlab("Offspring Acclimation Temperature (Degrees C)") + ylab("Deviance Residuals")

    Res.by.Mat = ggplot(data = data.frame("Mat" = as.factor(MMR.NMS.Data$Mother.Acc),
    "Res" = residuals(MMR.m3, type = "deviance")), aes(x = Mat, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "mediumseagreen") + my.theme.diag + 
    xlab("Dam Acclimation Temperature (Degrees C)") + ylab("Deviance Residuals")

    QQNorm.by.Off = ggplot(data.frame("Res" = residuals(MMR.m3, type = "deviance"), 
    "Fit" = predict(MMR.m3), "Off" = as.factor(MMR.NMS.Data$Offspring.Acc)), 
    aes(sample=Res, colour = Off)) + stat_qq() + 
    stat_qq_line() + my.theme + theme_bw() + 
    scale_colour_manual(values = c("mediumseagreen", "mediumorchid")) + theme(legend.position = "top") +
    guides(colour=guide_legend(title="Offspring Acclimation Temperature"))

    QQNorm.by.Mat = ggplot(data.frame("Res" = residuals(MMR.m3, type = "deviance"), 
    "Fit" = predict(MMR.m3), "Mat" = as.factor(MMR.NMS.Data$Mother.Acc)), 
    aes(sample=Res, colour = Mat)) + stat_qq() + 
    stat_qq_line() + my.theme + theme_bw() + 
    scale_colour_manual(values = c("royalblue", "gold2")) + theme(legend.position = "top") +
    guides(colour=guide_legend(title="Dam Acclimation Temperature"))

    grid.arrange(Res.by.Mass, Res.by.Off, Res.by.Mat, QQNorm.by.Off, QQNorm.by.Pat,
    nrow = 3, top = "Maximum Metabolic Rate Residuals") 

 ## Similar variance as above two models. Perhaps an interaction between mass and temperature categories
 ## could be considered? Calculating p-values.

    dwdata.MMR.NMS.m3 = data.frame(term = c(rownames(summary(MMR.m3)$coefficients)),
                       estimate = c(summary(MMR.m3)$coefficients[1:4]),
                       std.error = c(summary(MMR.m3)$coefficients[5:8]),
                       statistic = c(summary(MMR.m3)$coefficients[9:12]))

    p.values.m3= c()
    for (i in 1:nrow(dwdata.MMR.NMS.m3)){
      p.values.m3[i] = pt(abs(dwdata.MMR.NMS.m3$statistic[i]), df = 8, lower.tail = FALSE)
    } 
    dwdata.MMR.NMS.m3$p.values = round(p.values.m3, 3)

    dwdata.MMR.NMS.m3$term = c("(Intercept)", "log Mass", "Offspring Acclimation Temperature", 
    "Dam Acclimation Temperature")

## Comparing Model Outputs by Plots 

dwdata.MMR.NMS.m1 = dwdata.MMR.NMS.m1 %>% mutate(model = "Model 1")
dwdata.MMR.NMS.m2 = dwdata.MMR.NMS.m2 %>% mutate(model = "Model 2")
dwdata.MMR.NMS.m3 = dwdata.MMR.NMS.m3 %>% mutate(model = "Model 3")

MMR.NMS.Models = rbind(dwdata.MMR.NMS.m1, dwdata.MMR.NMS.m2, dwdata.MMR.NMS.m3)

MMR.NMS.Models$estimate = MMR.NMS.Models$estimate*8
MMR.NMS.Models$std.error = MMR.NMS.Models$std.error*8 

MMR.NMS.Models %>%
  dwplot(style = "distribution", show_intercept = FALSE) + theme_bw() + 
  theme(legend.justification = c(-0.1,-0.1), 
  legend.position = c(0,0.01)) + my.theme.dw + xlab("Coefficient") + ylab("") + 
  geom_vline(xintercept = 0, colour = "black", linetype = 2) + 
  scale_fill_manual(values = c("mediumseagreen", "mediumorchid", "cornflowerblue")) + 
  scale_colour_manual(values = c("black", "black", "black"))  + scale_x_continuous(limits = c(-8.0,8.0), 
  breaks = c(-8.0, -4.0, 0, 4.0, 8.0), label = c("-1.00", "-0.50", "0", "0.50", "1.00"))

ggsave("/home/joshk/git_repositories/BrookTrout_TG/Figures/log Maximum Metabolic Rate Model Comparisons.jpeg", dpi = 800, limitsize = TRUE)

MMR.NMS.Models$estimate = MMR.NMS.Models$estimate/8
MMR.NMS.Models$std.error = MMR.NMS.Models$std.error/8 

MMR.NMS.Fit.m1 = data.frame("Response" = c(log(MMR.NMS.Data$MMR)), "Fit" = c(predict(MMR.m1)), 
    "Mass" = c(log(MMR.NMS.Data$Mass)))
MMR.NMS.Fit.m2 = data.frame("Response" = c(log(MMR.NMS.Data$MMR)), "Fit" = c(predict(MMR.m2)), 
    "Mass" = c(log(MMR.NMS.Data$Mass)))
MMR.NMS.Fit.m3 = data.frame("Response" = c(log(MMR.NMS.Data$MMR)), "Fit" = c(predict(MMR.m3)), 
    "Mass" = c(log(MMR.NMS.Data$Mass)))

MMR.NMS.Fit.m1 = MMR.NMS.Fit.m1 %>% mutate(model = "Model 1")
MMR.NMS.Fit.m2 = MMR.NMS.Fit.m2 %>% mutate(model = "Model 2")
MMR.NMS.Fit.m3 = MMR.NMS.Fit.m3 %>% mutate(model = "Model 3")

MMR.NMS.Fit.All = rbind(MMR.NMS.Fit.m1, MMR.NMS.Fit.m2, MMR.NMS.Fit.m3)
MMR.NMS.Fits = ggplot(data = MMR.NMS.Fit.All, aes(x = Fit, y = Response, fill = model, colour = model,
  linetype = model)) + theme_bw() + my.theme + geom_point(pch=21, size = 3, alpha = 0.5) + 
  geom_smooth(method = "lm") + xlab("Predicted log Resting Metabolic Rate \n (log mg O2/h)") + 
  ylab("Measured log Resting \n Metabolic Rate \n (log mg O2/h)") + 
  scale_fill_manual(values = c("mediumseagreen", 
  "mediumorchid", "cornflowerblue")) + scale_colour_manual(values = c("black", 
  "black", "black")) + theme(legend.title=element_blank())

print(MMR.NMS.Fits)

ggsave("/home/joshk/git_repositories/BrookTrout_TG/Figures/log Maximum Metabolic Rate Model Fit Comparisons.jpeg", dpi = 800, limitsize = TRUE)

# Very subtle differences between models.

###Testing Model Assumptions for Absolute Aerobic Scope: Model 1.

    AAS.m1 = lmer(AAS ~ log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID), 
        REML = FALSE, data = AAS.NMS.Data)

    summary(AAS.m1)$varcor
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.0000  
##  Father.ID (Intercept) 0.0000  
##  Residual              0.2058
# Maternal ID explaining zero variance.

    Res.Alone = ggplot(data = data.frame("Res" = residuals(AAS.m1, type = "deviance")), 
    aes(x = 1:nrow(AAS.NMS.Data), y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
    ylab("Deviance Residuals")

    Res.by.Fit = ggplot(data = data.frame("Res" = residuals(AAS.m1, type = "deviance"), 
    "Fit" = predict(AAS.m1)), aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
    ylab("Deviance Residuals")

    Res.by.Resp = ggplot(data = data.frame("Res" = residuals(AAS.m1, type = "deviance"), 
    "Resp" = AAS.NMS.Data$AAS), aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
    xlab("Resting Metabolic Rate") + ylab("Deviance Residuals")

    QQNorm = ggplot(data.frame("Res" = residuals(AAS.m1, type = "deviance"), "Fit" = predict(AAS.m1)),
    aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = 
    "Absolute Aerobic Scope Residuals")  

    # Residuals normal, fan in residuals vs fitteds.

    Resid.by.Mass = ggplot(data = data.frame("Res" = residuals(AAS.m1, type = "deviance"), 
    "Mass" = log(AAS.NMS.Data$Mass)), aes(x = Mass, y = Res)) + my.theme + theme_bw() + geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + 
    xlab("log Mass (g)") + ylab("Deviance Residuals")

    print(Resid.by.Mass)

    # Weight by inverse of mass - conducted below and compared to original model.

    AAS.m1b = lmer(AAS ~ log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID), 
        REML = FALSE, weights = 1/log(Mass), data = AAS.NMS.Data)

R2vals = c(unlist(lapply(list(AAS.m1, AAS.m1b), r.squaredGLMM)))

Model.Comparison = data.frame("Weighted" = c("No", "Yes"), "Marginal R2" = c(R2vals[c(1,3)]),
    "Conditional R2" = c(R2vals[c(2,4)]), "AIC" = c(summary(AAS.m1)$AICtab[1], summary(AAS.m1b)$AICtab[1]))

print(Model.Comparison)
##   Weighted Marginal.R2 Conditional.R2      AIC
## 1       No   0.4030555      0.4030555 -42.7185
## 2      Yes   0.4419980      0.4419980 -52.5377
    # Both variance explained and AIC improved. Assessing mass by residual structure

    Mass.by.Res = ggplot(data = data.frame("Res" = residuals(AAS.m1b, type = "pearson"),
    "Mass" = log(AAS.NMS.Data$Mass)), 
    aes(x = Mass, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + 
    xlab("log Mass (g)") + ylab("Deviance Residuals")

    print(Mass.by.Res)

###Re-running Model Iterations for Absolute Aerobic Scope: Application of 1/log Mass Weighting.

    AAS.model = lmer(AAS ~ Mother.Acc*Father.Acc*Offspring.Acc + log(Mass) + (1|Mother.ID) +(1|Father.ID), 
    REML = FALSE, weights = 1/log(Mass), data = AAS.NMS.Data)   

    AIC.Results.AAS = subset(dredge(AAS.model, rank = "AICc", 
        extra = list(R2 = function(x) r.squaredGLMM(x))), delta < 5)

Select.Models.AAS = c()
Delta.AAS = c()
Response = c()
Marg.R2 = c() 
Cond.R2 = c()

  for (i in 1:length(attr(AIC.Results.AAS, "model.calls"))){
    Select.Models.AAS[i] = gsub(".*~ ", "", paste(attr(AIC.Results.AAS, "model.calls")[[i]])[2])
    Delta.AAS[i] = AIC.Results.AAS$delta[i]
    Response[i] = gsub(" ~.*", "", paste(attr(AIC.Results.AAS, "model.calls")[[i]])[2])
    Marg.R2[i] = AIC.Results.AAS$R21[i]
    Cond.R2[i] = AIC.Results.AAS$R22[i]
    if (i == length(attr(AIC.Results.AAS, "model.calls"))){
      Table = data.frame("Reponse.Type" = "Mass Independant", "Response" = Response, 
      "Delta.AIC" = Delta.AAS, "Marginal.R2" = Marg.R2, "Conditional.R2" = Cond.R2,
      "Equation" = Select.Models.AAS)
      assign(paste("AIC.Selection", "AAS", sep = "_"), Table)
    }
  }

    AIC.Selection.NMS = rbind(AIC.Selection_NMS_1, AIC.Selection_MMR, AIC.Selection_AAS, AIC.Selection_NMS_4)

    AIC.Selection = rbind(AIC.Selection.MS, AIC.Selection.NMS)
    print(AIC.Selection, right = F)
##    Reponse.Type     Response    Delta.AIC Marginal.R2 Conditional.R2
## 1  Mass-Specific    log(RMR.MS) 0.0000000 0.09689341  0.23258538    
## 2  Mass-Specific    log(RMR.MS) 0.4320103 0.05317331  0.24043771    
## 3  Mass-Specific    log(RMR.MS) 1.2633111 0.11949983  0.22942058    
## 4  Mass-Specific    log(RMR.MS) 1.6197476 0.07652765  0.23332101    
## 5  Mass-Specific    log(RMR.MS) 2.1592505 0.09688383  0.23257954    
## 6  Mass-Specific    log(RMR.MS) 3.4158810 0.11948859  0.22852419    
## 7  Mass-Specific    log(RMR.MS) 3.4339163 0.11966268  0.23078182    
## 8  Mass-Specific    log(RMR.MS) 3.4465742 0.11949854  0.22941555    
## 9  Mass-Specific    log(RMR.MS) 3.7625995 0.07663967  0.23491827    
## 10 Mass-Specific    log(MMR.MS) 0.0000000 0.12081940  0.14201563    
## 11 Mass-Specific    log(MMR.MS) 1.1337870 0.12482063  0.14707607    
## 12 Mass-Specific    log(MMR.MS) 1.6198084 0.07304761  0.14942805    
## 13 Mass-Specific    log(MMR.MS) 1.6632302 0.12283706  0.14581838    
## 14 Mass-Specific    log(MMR.MS) 2.7260488 0.07701051  0.15563521    
## 15 Mass-Specific    log(MMR.MS) 2.7981663 0.12640953  0.14714116    
## 16 Mass-Specific    log(MMR.MS) 2.8091275 0.12687437  0.15094298    
## 17 Mass-Specific    log(MMR.MS) 3.1451250 0.12529971  0.14712299    
## 18 Mass-Specific    log(MMR.MS) 4.4679396 0.07829618  0.15562601    
## 19 Mass-Specific    log(MMR.MS) 4.5445400 0.12823397  0.15073801    
## 20 Mass-Specific    log(MMR.MS) 4.7897724 0.12703622  0.14722626    
## 21 Mass-Specific    log(MMR.MS) 4.8318359 0.12737623  0.15100909    
## 22 Mass-Specific    log(AAS.MS) 0.0000000 0.06030224  0.06030224    
## 23 Mass-Specific    log(AAS.MS) 1.6058206 0.06331924  0.06331924    
## 24 Mass-Specific    log(AAS.MS) 1.6547206 0.06304827  0.06304827    
## 25 Mass-Specific    log(AAS.MS) 2.7150419 0.06921105  0.06921105    
## 26 Mass-Specific    log(AAS.MS) 3.2960644 0.06600676  0.06600676    
## 27 Mass-Specific    log(AAS.MS) 3.8171171 0.06312377  0.06312377    
## 28 Mass-Specific    log(AAS.MS) 4.4085037 0.07201229  0.07201229    
## 29 Mass Independant log(RMR)    0.0000000 0.13398490  0.18729777    
## 30 Mass Independant log(RMR)    1.2602254 0.14177631  0.20537453    
## 31 Mass Independant log(RMR)    1.6073708 0.15042976  0.21764843    
## 32 Mass Independant log(RMR)    1.6474074 0.12869909  0.17308051    
## 33 Mass Independant log(RMR)    2.7321064 0.13274864  0.17071098    
## 34 Mass Independant log(RMR)    2.8878776 0.13578794  0.19007496    
## 35 Mass Independant log(RMR)    3.2743628 0.14463523  0.20289714    
## 36 Mass Independant log(RMR)    4.0588194 0.14000948  0.18880081    
## 37 Mass Independant log(RMR)    4.0950899 0.09988873  0.21612058    
## 38 Mass Independant log(RMR)    4.5640317 0.14806374  0.20017371    
## 39 Mass Independant log(RMR)    4.8488012 0.13927880  0.19821267    
## 40 Mass Independant log(MMR)    0.0000000 0.51073885  0.52414430    
## 41 Mass Independant log(MMR)    1.1519770 0.50259058  0.51232882    
## 42 Mass Independant log(MMR)    1.6571485 0.51244418  0.52665944    
## 43 Mass Independant log(MMR)    2.3452891 0.51490234  0.52733792    
## 44 Mass Independant log(MMR)    2.7795841 0.50375882  0.51398621    
## 45 Mass Independant log(MMR)    3.2968819 0.50262005  0.51249178    
## 46 Mass Independant log(MMR)    3.4572491 0.50709834  0.51576702    
## 47 Mass Independant log(MMR)    4.6305941 0.48890967  0.49534628    
## 48 Mass Independant log(MMR)    4.8358232 0.50402584  0.51443503    
## 49 Mass Independant log(MMR)    4.9419345 0.50379582  0.51418385    
## 50 Mass Independant AAS         0.0000000 0.44199801  0.44199801    
## 51 Mass Independant AAS         2.1304077 0.44178239  0.44178239    
## 52 Mass Independant AAS         2.1610466 0.44220786  0.44220786    
## 53 Mass Independant AAS         2.7171750 0.44576430  0.44576430    
## 54 Mass Independant AAS         4.3196162 0.44198897  0.44198897    
## 55 Mass Independant AAS         4.3650932 0.44224301  0.44224301    
## 56 Mass Independant AAS         4.9188507 0.44606639  0.44606639    
## 57 Mass Independant CTM         0.0000000 0.46835156  0.50879302    
##    Equation                                                                                                                                  
## 1  Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 2  Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                         
## 3  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                               
## 4  Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 5  Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                 
## 6  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                       
## 7  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                    
## 8  Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                    
## 9  Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                 
## 10 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 11 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                               
## 12 Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                         
## 13 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                 
## 14 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 15 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                    
## 16 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                    
## 17 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                       
## 18 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                 
## 19 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc + Mother.Acc:Offspring.Acc         
## 20 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc + Mother.Acc:Offspring.Acc            
## 21 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc + Father.Acc:Offspring.Acc            
## 22 Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                         
## 23 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 24 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                            
## 25 Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                                 
## 26 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                               
## 27 Father.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                                 
## 28 Father.Acc + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                    
## 29 log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                                
## 30 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 31 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                     
## 32 Father.Acc + log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                   
## 33 Father.Acc + log(Mass) + Mother.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                                           
## 34 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                   
## 35 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                        
## 36 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                           
## 37 log(Mass) + (1 | Mother.ID) + (1 | Father.ID)                                                                                             
## 38 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc + Mother.Acc:Offspring.Acc
## 39 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                        
## 40 log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                             
## 41 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 42 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 43 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                     
## 44 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                   
## 45 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                     
## 46 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                        
## 47 log(Mass) + (1 | Mother.ID) + (1 | Father.ID)                                                                                             
## 48 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Mother.Acc                           
## 49 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                        
## 50 log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                             
## 51 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 52 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                                
## 53 log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                                     
## 54 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID)                                                   
## 55 Father.Acc + log(Mass) + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Father.Acc:Offspring.Acc                                     
## 56 Father.Acc + log(Mass) + Mother.Acc + Offspring.Acc + (1 | Mother.ID) + (1 | Father.ID) + Mother.Acc:Offspring.Acc                        
## 57 Offspring.Acc
    write.csv(AIC.Selection, "Model_Outputs_AICc.csv", row.names = FALSE)

###Reassessing Model Assumptions for Absolute Aerobic Scope: Model 1.

    AAS.m1 = lmer(AAS ~ log(Mass) + Offspring.Acc +
    (1|Mother.ID) + (1|Father.ID), REML = FALSE,  weights = 1/log(Mass), 
    data = AAS.NMS.Data)

    summary(AAS.m1)$varcor
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.00000 
##  Father.ID (Intercept) 0.00000 
##  Residual              0.19642
    # Neither paternal nor maternal ID explain residual variance. 

    Res.Alone = ggplot(data = data.frame("Res" = residuals(AAS.m1, type = "pearson")), 
    aes(x = 1:nrow(AAS.NMS.Data), y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
    ylab("Pearson Residuals")

    Res.by.Fit = ggplot(data = data.frame("Res" = residuals(AAS.m1, type = "pearson"), "Fit" = predict(AAS.m1)), 
    aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
    ylab("Pearson Residuals")

    Res.by.Resp = ggplot(data = data.frame("Res" = residuals(AAS.m1, type = "pearson"), "Resp" = 
    AAS.NMS.Data$AAS), aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
    xlab("Resting Metabolic Rate") + ylab("Pearson Residuals")

    QQNorm = ggplot(data.frame("Res" = residuals(AAS.m1, type = "pearson"), "Fit" = predict(AAS.m1)),
    aes(sample=Res)) + stat_qq(colour = "gold") + stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = 
    "Absolute Aerobic Scope Residuals")  

    # Well behaved.

    qqPlot(residuals(AAS.m1, type = "pearson"), ylab = "Pearson Residuals")

## 102 206 
##  76 152
    # Normal pearson residuals.

    # Assessing residuals by predictors.
 
    Res.by.Mass = ggplot(data = data.frame("Res" = residuals(AAS.m1, type = "pearson"), 
    "Mass" = log(AAS.NMS.Data$Mass)), aes(x = Mass, y = Res)) + theme_bw() + my.theme.diag + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + 
    xlab("log Mass (g)") + ylab("Deviance Residuals")

    Res.by.Off = ggplot(data = data.frame("Off" = as.factor(AAS.NMS.Data$Offspring.Acc),
    "Res" = residuals(AAS.m1, type = "pearson")), aes(x = Off, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "royalblue") + my.theme.diag + 
    xlab("Offspring Acclimation Temperature (Degrees C)") + ylab("Pearson Residuals")

    grid.arrange(Res.by.Mass, Res.by.Off, nrow = 2, top = 
    "Absolute Aerobic Scope Residuals")  

    # Very subtle heterogeneity, but not sufficient for adjusted variance structure. 
    # Calculating p-values.

    dwdata.AAS.NMS.m1 = data.frame(term = c(rownames(summary(AAS.m1)$coefficients)),
                       estimate = c(summary(AAS.m1)$coefficients[1:3]),
                       std.error = c(summary(AAS.m1)$coefficients[4:6]),
                       statistic = c(summary(AAS.m1)$coefficients[7:9]))

    p.values.m1= c()
    for (i in 1:nrow(dwdata.AAS.NMS.m1)){
     p.values.m1[i] = pt(abs(dwdata.AAS.NMS.m1$statistic[i]), df = 8, lower.tail = FALSE)
    } 
    dwdata.AAS.NMS.m1$p.values = round(p.values.m1, 3)

    dwdata.AAS.NMS.m1$term = c("(Intercept)", "log Mass", 
    "Offspring Acclimation Temperature" 
    )

###Testing Model Assumptions for Absolute Aerobic Scope: Model 2.

    AAS.m2 = lmer(AAS ~ log(Mass) + Offspring.Acc + Mother.Acc + 
    (1|Mother.ID) + (1|Father.ID), REML = FALSE,  weights = 1/log(Mass), 
    data = AAS.NMS.Data)

    summary(AAS.m2)$varcor
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.00000 
##  Father.ID (Intercept) 0.00000 
##  Residual              0.19639
    # Neither paternal nor maternal ID explain residual variance, as in the previous model. 

    Res.Alone = ggplot(data = data.frame("Res" = residuals(AAS.m2, type = "pearson")), 
    aes(x = 1:nrow(AAS.NMS.Data), y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
    ylab("Pearson Residuals")

    Res.by.Fit = ggplot(data = data.frame("Res" = residuals(AAS.m2, type = "pearson"), 
    "Fit" = predict(AAS.m2)), aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
    ylab("Pearson Residuals")

    Res.by.Resp = ggplot(data = data.frame("Res" = residuals(AAS.m2, type = "pearson"), 
    "Resp" = AAS.NMS.Data$AAS), aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
    xlab("Resting Metabolic Rate") + ylab("Pearson Residuals")

    QQNorm = ggplot(data.frame("Res" = residuals(AAS.m2, type = "pearson"), 
    "Fit" = predict(AAS.m2)), aes(sample=Res)) + stat_qq(colour = "gold") + 
    stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = 
    "Absolute Aerobic Scope Residuals")  

    # Similarly well behaved.

    qqPlot(residuals(AAS.m2, type = "pearson"), ylab = "Pearson Residuals")

## 206 102 
## 152  76
    # Normal pearson residuals, though with slight tails.

    # Assessing residuals by predictors.
 
    Res.by.Mass = ggplot(data = data.frame("Res" = residuals(AAS.m2, type = "pearson"), 
    "Mass" = log(AAS.NMS.Data$Mass)), aes(x = Mass, y = Res)) + theme_bw() + my.theme.diag + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + 
    xlab("log Mass (g)") + ylab("Deviance Residuals")

    Res.by.Off = ggplot(data = data.frame("Off" = as.factor(AAS.NMS.Data$Offspring.Acc),
    "Res" = residuals(AAS.m2, type = "pearson")), aes(x = Off, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "gold2") + my.theme.diag + 
    xlab("Offspring Acclimation Temperature (Degrees C)") + ylab("Pearson Residuals")

    Res.by.Dam = ggplot(data = data.frame("Dam" = as.factor(AAS.NMS.Data$Mother.Acc),
    "Res" = residuals(AAS.m2, type = "pearson")), aes(x = Dam, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "mediumorchid") + my.theme.diag + 
    xlab("Dam Acclimation Temperature (Degrees C)") + ylab("Pearson Residuals")

    layout = rbind(c(1,1,2,2), c(NA,3,3,NA))
    grid.arrange(Res.by.Mass, Res.by.Off, Res.by.Dam, nrow = 2, top = 
    "Absolute Aerobic Scope Residuals", layout_matrix = layout)  

    # Very subtle heteroskedasticity across dam acclimation temperature, 
    # but again, likely not sufficient to change variance structure. 

    # Calculating p-values.

    dwdata.AAS.NMS.m2 = data.frame(term = c(rownames(summary(AAS.m2)$coefficients)),
                       estimate = c(summary(AAS.m2)$coefficients[1:4]),
                       std.error = c(summary(AAS.m2)$coefficients[5:8]),
                       statistic = c(summary(AAS.m2)$coefficients[9:12]))

    p.values.m2 = c()
    for (i in 1:nrow(dwdata.AAS.NMS.m2)){
     p.values.m2[i] = pt(abs(dwdata.AAS.NMS.m2$statistic[i]), df = 8, lower.tail = FALSE)
    } 
    dwdata.AAS.NMS.m2$p.values = round(p.values.m2, 3)

    dwdata.AAS.NMS.m2$term = c("(Intercept)", "log Mass", 
    "Offspring Acclimation Temperature", "Dam Acclimation Temperature" 
    )

###Testing Model Assumptions for Absolute Aerobic Scope: Model 3.

    AAS.m3 = lmer(AAS ~ log(Mass) + Offspring.Acc + Father.Acc + 
    (1|Mother.ID) + (1|Father.ID), REML = FALSE,  weights = 1/log(Mass), 
    data = AAS.NMS.Data)

    summary(AAS.m3)$varcor
##  Groups    Name        Std.Dev.
##  Mother.ID (Intercept) 0.00000 
##  Father.ID (Intercept) 0.00000 
##  Residual              0.19641
    # Again, neither paternal nor maternal ID explain residual variance. 

    Res.Alone = ggplot(data = data.frame("Res" = residuals(AAS.m3, type = "pearson")), 
    aes(x = 1:nrow(AAS.NMS.Data), y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) + xlab("Row Number") + 
    ylab("Pearson Residuals")

    Res.by.Fit = ggplot(data = data.frame("Res" = residuals(AAS.m3, type = "pearson"), 
    "Fit" = predict(AAS.m3)), aes(x = Fit, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + xlab("Fitted Values") + 
    ylab("Pearson Residuals")

    Res.by.Resp = ggplot(data = data.frame("Res" = residuals(AAS.m3, type = "pearson"), 
    "Resp" = AAS.NMS.Data$AAS), aes(x = Resp, y = Res)) + my.theme + theme_bw() + 
    geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) + 
    xlab("Resting Metabolic Rate") + ylab("Pearson Residuals")

    QQNorm = ggplot(data.frame("Res" = residuals(AAS.m3, type = "pearson"), 
    "Fit" = predict(AAS.m3)), aes(sample=Res)) + stat_qq(colour = "gold") + 
    stat_qq_line() + my.theme + theme_bw()

    grid.arrange(Res.Alone, Res.by.Fit, Res.by.Resp, QQNorm, nrow = 2, top = 
    "Absolute Aerobic Scope Residuals")  

    # No distinct residual patterns.

    qqPlot(residuals(AAS.m3, type = "pearson"), ylab = "Pearson Residuals")

## 102 206 
##  76 152
    # Normal pearson residuals, as before.

    # Again, assessing residuals by predictors.
 
    Res.by.Mass = ggplot(data = data.frame("Res" = residuals(AAS.m3, type = "pearson"), 
    "Mass" = log(AAS.NMS.Data$Mass)), aes(x = Mass, y = Res)) + theme_bw() + my.theme.diag + 
    geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) + 
    xlab("log Mass (g)") + ylab("Deviance Residuals")

    Res.by.Off = ggplot(data = data.frame("Off" = as.factor(AAS.NMS.Data$Offspring.Acc),
    "Res" = residuals(AAS.m3, type = "pearson")), aes(x = Off, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "gold2") + my.theme.diag + 
    xlab("Offspring Acclimation Temperature (Degrees C)") + ylab("Pearson Residuals")

    Res.by.Pat = ggplot(data = data.frame("Pat" = as.factor(AAS.NMS.Data$Mother.Acc),
    "Res" = residuals(AAS.m3, type = "pearson")), aes(x = Pat, y = Res)) + 
    theme_bw() + geom_boxplot(fill = "mediumorchid") + my.theme.diag + 
    xlab("Dam Acclimation Temperature (Degrees C)") + ylab("Pearson Residuals")

    layout = rbind(c(1,1,2,2), c(NA,3,3,NA))
    grid.arrange(Res.by.Mass, Res.by.Off, Res.by.Pat, nrow = 2, top = 
    "Absolute Aerobic Scope Residuals", layout_matrix = layout)  

    # Again, very subtle heteroskedasticity across both sire and offspring acclimation 
    # temperature, though too minor to warrant concern.

    # Calculating p-values.

    dwdata.AAS.NMS.m3 = data.frame(term = c(rownames(summary(AAS.m3)$coefficients)),
                       estimate = c(summary(AAS.m2)$coefficients[1:4]),
                       std.error = c(summary(AAS.m2)$coefficients[5:8]),
                       statistic = c(summary(AAS.m2)$coefficients[9:12]))

    p.values.m3 = c()
    for (i in 1:nrow(dwdata.AAS.NMS.m3)){
     p.values.m3[i] = pt(abs(dwdata.AAS.NMS.m3$statistic[i]), df = 8, lower.tail = FALSE)
    } 
    dwdata.AAS.NMS.m3$p.values = round(p.values.m3, 3)

    dwdata.AAS.NMS.m3$term = c("(Intercept)", "log Mass", 
    "Offspring Acclimation Temperature", "Sire Acclimation Temperature" 
    )

###Comparing Absolute Aerobic Scope Models by Plot.

dwdata.AAS.NMS.m1 = dwdata.AAS.NMS.m1 %>% mutate(model = "Model 1")
dwdata.AAS.NMS.m2 = dwdata.AAS.NMS.m2 %>% mutate(model = "Model 2")
dwdata.AAS.NMS.m3 = dwdata.AAS.NMS.m3 %>% mutate(model = "Model 3")

AAS.NMS.Models = rbind(dwdata.AAS.NMS.m1, dwdata.AAS.NMS.m2, dwdata.AAS.NMS.m3)

AAS.NMS.Models$estimate = AAS.NMS.Models$estimate*10
AAS.NMS.Models$std.error = AAS.NMS.Models$std.error*10 

AAS.NMS.Models %>%
  dwplot(style = "distribution", show_intercept = FALSE) + theme_bw() + theme(legend.justification = c(-0.1,-0.1), 
  legend.position = c(0,0.01)) + my.theme.dw + xlab("Coefficient") + ylab("") + geom_vline(xintercept = 0, 
  colour = "black", linetype = 2) + scale_fill_manual(values = c("mediumseagreen", "mediumorchid", "cornflowerblue")) + 
  scale_colour_manual(values = c("black", "black", "black"))  + scale_x_continuous(limits = c(-10.0,10.0), 
  breaks = c(-10.0, -5.0, 0, 5.0, 10.0), label = c("-1.00", "-0.50", "0", "0.50", "1.00"))

ggsave("/home/joshk/git_repositories/BrookTrout_TG/Figures/Absolute Aerobic Scope Model Comparisons.jpeg", dpi = 800, limitsize = TRUE)

AAS.NMS.Models$estimate = AAS.NMS.Models$estimate/10
AAS.NMS.Models$std.error = AAS.NMS.Models$std.error/10 

AAS.NMS.Fit.m1 = data.frame("Response" = c(AAS.NMS.Data$AAS), "Fit" = c(predict(AAS.m1)), 
    "Mass" = c(log(AAS.NMS.Data$Mass)))
AAS.NMS.Fit.m2 = data.frame("Response" = c(AAS.NMS.Data$AAS), "Fit" = c(predict(AAS.m2)), 
    "Mass" = c(log(AAS.NMS.Data$Mass)))
AAS.NMS.Fit.m3 = data.frame("Response" = c(AAS.NMS.Data$AAS), "Fit" = c(predict(AAS.m3)), 
    "Mass" = c(log(AAS.NMS.Data$Mass)))

AAS.NMS.Fit.m1 = AAS.NMS.Fit.m1 %>% mutate(model = "Model 1")
AAS.NMS.Fit.m2 = AAS.NMS.Fit.m2 %>% mutate(model = "Model 2")
AAS.NMS.Fit.m3 = AAS.NMS.Fit.m3 %>% mutate(model = "Model 3")

AAS.NMS.Fit.All = rbind(AAS.NMS.Fit.m1, AAS.NMS.Fit.m2, AAS.NMS.Fit.m3)
AAS.NMS.Fits = ggplot(data = AAS.NMS.Fit.All, aes(x = Fit, y = Response, fill = model, colour = model,
  linetype = model)) + theme_bw() + my.theme + geom_point(pch=21, size = 3, alpha = 0.5) + 
  geom_smooth(method = "lm") + xlab("Predicted Absolute \n Aerobic Scope (mg O2/h)") + 
  ylab("Measured Absolute \n Aerobic Scope (mg O2/h)") + scale_fill_manual(values = c("mediumseagreen", 
  "mediumorchid", "cornflowerblue")) + scale_colour_manual(values = c("black", 
  "black", "black")) + theme(legend.title=element_blank())

print(AAS.NMS.Fits)

ggsave("/home/joshk/git_repositories/BrookTrout_TG/Figures/Absolute Aerobic Scope Model Fit Comparisons.jpeg", dpi = 800, limitsize = TRUE)

###Amalgamating Mass Independant Results.

require('plyr')

R2.Dataframe = ldply(lapply(list(RMR.m1, RMR.m2, RMR.m3, 
        MMR.m1, MMR.m2, MMR.m3, 
        AAS.m1, AAS.m2, AAS.m3),
        r.squaredGLMM), data.frame)

RMR.NMS.DF1 = subset(RMR.NMS.Models, model == "Model 1") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 1")%>% 
  mutate("Marginal R2" = R2.Dataframe[1,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[1,2])%>%
  mutate("AIC" = summary(RMR.m1)$AICtab[1])
RMR.NMS.DF2 = subset(RMR.NMS.Models, model == "Model 2") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 2")%>% 
  mutate("Marginal R2" = R2.Dataframe[2,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[2,2]) %>%
  mutate("AIC" = summary(RMR.m2)$AICtab[1]) 
RMR.NMS.DF3 = subset(RMR.NMS.Models, model == "Model 3") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 3")%>% 
  mutate("Marginal R2" = R2.Dataframe[3,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[3,2])%>% 
  mutate("AIC" = summary(RMR.m3)$AICtab[1])

RMR.NMS.DF = rbind(RMR.NMS.DF1, RMR.NMS.DF2, RMR.NMS.DF3)

RMR.NMS.DF$Response = "Resting Metabolic Rate" 
RMR.NMS.DF$Response = "Resting Metabolic Rate" 
RMR.NMS.DF$Response = "Resting Metabolic Rate" 

MMR.NMS.DF1 = subset(MMR.NMS.Models, model == "Model 1") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 1")%>% 
  mutate("Marginal R2" = R2.Dataframe[4,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[4,2]) %>% 
  mutate("AIC" = summary(MMR.m1)$AICtab[1])
MMR.NMS.DF2 = subset(MMR.NMS.Models, model == "Model 2") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 2")%>% 
  mutate("Marginal R2" = R2.Dataframe[5,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[5,2]) %>% 
  mutate("AIC" = summary(MMR.m2)$AICtab[1])
MMR.NMS.DF3 = subset(MMR.NMS.Models, model == "Model 3") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 3")%>% 
  mutate("Marginal R2" = R2.Dataframe[6,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[6,2]) %>% 
  mutate("AIC" = summary(MMR.m3)$AICtab[1])

MMR.NMS.DF = rbind(MMR.NMS.DF1, MMR.NMS.DF2, MMR.NMS.DF3)

MMR.NMS.DF$Response = "Maximum Metabolic Rate" 
MMR.NMS.DF$Response = "Maximum Metabolic Rate" 
MMR.NMS.DF$Response = "Maximum Metabolic Rate" 

AAS.NMS.DF1 = subset(AAS.NMS.Models, model == "Model 1") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 1")%>% 
  mutate("Marginal R2" = R2.Dataframe[7,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[7,2]) %>% 
  mutate("AIC" = summary(AAS.m1)$AICtab[1])

AAS.NMS.DF2 = subset(AAS.NMS.Models, model == "Model 2") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 2")%>% 
  mutate("Marginal R2" = R2.Dataframe[8,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[8,2]) %>% 
  mutate("AIC" = summary(AAS.m2)$AICtab[1])

AAS.NMS.DF3 = subset(AAS.NMS.Models, model == "Model 3") %>% mutate("df" = 8)%>%
  mutate("Model Number" = "Model 3")%>% 
  mutate("Marginal R2" = R2.Dataframe[9,1])%>% 
  mutate("Conditional R2" = R2.Dataframe[9,2]) %>% 
  mutate("AIC" = summary(AAS.m3)$AICtab[1])

AAS.NMS.DF = rbind(AAS.NMS.DF1, AAS.NMS.DF2, AAS.NMS.DF3)

AAS.NMS.DF$Response = "Absolute Aerobic Scope" 
AAS.NMS.DF$Response = "Absolute Aerobic Scope" 
AAS.NMS.DF$Response = "Absolute Aerobic Scope" 

AllData.NMS = rbind(RMR.NMS.DF, MMR.NMS.DF, AAS.NMS.DF) 

AllData.NMS = AllData.NMS %>% select("Response", "Model Number", "term", "estimate", "std.error", "statistic",
                                     "df", "p.values", "AIC", 
                             "Marginal R2", "Conditional R2")
colnames(AllData.NMS) = c("Response Variable", "Model Number", "Model Parameter", "Estimate (beta)",
                          "Std.Error", "T Statistic", "DF", "p", 
                      "AIC", "Marginal R2", "Conditional R2")

New.Parameters = c()
for(i in 1:nrow(AllData.NMS)){
    New.Parameters[i] = gsub("\n", "", AllData.NMS$"Model Parameter"[i])
    }

AllData.NMS$"Model Parameter" = New.Parameters
write.csv(AllData.NMS, "Mass Independant Model Output Results.csv")